{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 任务三、EM算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 前言\n",
    "EM算法是机器学习十大算法之一，它很简单，但是也同样很有深度，简单是因为它就分两步求解问题，\n",
    "- E步：求期望（expectation）\n",
    "- M步：求极大（maximization)\n",
    "\n",
    "深度在于它的数学推理涉及到比较繁杂的概率公式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EM算法引入\n",
    "概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是<font color=\"#F00\" size = \"4px\">当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，</font>我们讨论极大似然估计，极大后验概率估计与其类似。\n",
    "参考统计学习方法书中的一个例子来引入EM算法， 假设有3枚硬币，分别记做A、B、C，这些硬币正面出现的概率分别是$\\pi$、$p$、$q$，进行如下实验：\n",
    "- 先掷硬币A，根据结果选出硬币B和硬币C，正面选硬币B，反面选硬币C\n",
    "- 通过选择出的硬币，掷硬币的结果出现正面为1，反面为0\n",
    "如此独立地重复n次实验，我们当前规定n=10，则10次的结果如下所示：\n",
    "$$\n",
    "1,1,0,1,0,0,1,0,1,1\n",
    "$$\n",
    "假设只通过观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三个硬币出现正面的概率？\n",
    "我们来构建这样一个三硬币模型：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "P(y|\\theta) &=\\sum_{z}P(y,z|\\theta)=\\sum_{z}P(z|\\theta)P(y|z,\\theta) \\\\\n",
    "  &=\\pi p^{y}(1-p)^{1-y}+(1-\\pi)q^{y}(1-q)^{1-y}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "- 若$y=1$，表示这此看到的是正面，这个正面有可能是B的正面，也可能是C的正面，则$P(1|\\theta)=\\pi p+(1-\\pi)q$\n",
    "- 若$y=0$，则$P(0|\\theta)=\\pi (1-p)+(1-\\pi)(1-q)$\n",
    "\n",
    "y是观测变量，表示一次观测结果是1或0，z是隐藏变量，表示掷硬币A的结果，这个是观测不到结果的，$\\theta=(\\pi,p,q)$表示模型参数，将观测数据表示为$Y=(Y_1,Y_2,...,Y_n)^{T}$，未观测的数据表示为$Z=(Z_1,Z_2,...,Z_n)^{T}$，则观测函数的似然函数是：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "P(Y|\\theta)&=\\sum_{Z}P(Z|\\theta)P(Y|Z,\\theta)\\\\\n",
    "&=\\prod_{i=0} ( \\pi p^{y_i}(1-p)^{1-y_{i}}+(1-\\pi)q^{y_{i}}(1-q)^{1-y_{i}})\n",
    "\\end{aligned}\n",
    "$$\n",
    "考虑求模型参数$\\theta=(\\pi,p,q)$的极大似然估计，即：\n",
    "$$\n",
    "\\hat{\\theta}=arg\\max_{\\theta}logP(Y|\\theta)\n",
    "$$\n",
    "这个问题没有解析解，只有通过迭代方法来求解，EM算法就是可以用于求解这个问题的一种迭代算法，下面给出EM算法的迭代过程：\n",
    "- 首先选取初始值，记做$\\theta^{0}=(\\pi^{0},p^{0},q^{0})$，第i次的迭代参数的估计值为$\\theta^{i}=(\\pi^{i},p^{i},q^{i})$\n",
    "- E步：计算在模型参数$\\pi^{i}，p^{i}，q^{i}$下观测变量$y_i$来源于硬币B的概率：\n",
    "$$\n",
    "\\mu^{i+1}=\\frac{\\pi^{i}(p^{i})^{y_i}(1-p^i)^{1-y_i}}{\\pi^{i}(p^{i})^{y_i}(1-p^i)^{1-y_i}+(1-\\pi^{i})(q^{i})^{y_i}(1-p^i)^{1-y_i}}\n",
    "$$\n",
    "备注一下：这个公式的分母是$P(Y|\\theta)$，分子表示是来源与B硬币的概率。\n",
    "\n",
    "- M步：计算模型参数的新估计值：\n",
    "$$\n",
    "\\pi^{i+1}=\\frac{1}{n}\\sum_{j=1}^{n}\\mu_{j}^{i+1} \n",
    "$$\n",
    "因为B硬币A硬币出现正面的结果，所以A硬币概率就是$\\mu_{j}$的平均值。\n",
    "$$\n",
    "p^{i+1}=\\frac{\\sum_{j=1}^{n}\\mu_{j}^{i+1}y_j}{\\sum_{j=1}^{n}\\mu_{j}^{i+1}}\n",
    "$$\n",
    "分子乘以$y_{i}$，所以其实是计算B硬币出现正面的概率。\n",
    "$$\n",
    "q^{i+1}=\\frac{\\sum_{j=1}^{n}(1-\\mu_{j}^{i+1})y_j}{\\sum_{j=1}^{n}(1-\\mu_{j}^{i+1})}\n",
    "$$\n",
    "$(1-\\mu_{j}^{i+1})$表示出现C硬币的概率。\n",
    "\n",
    "闭环形成，从$P(Y|\\theta)$ 到 $\\pi、p、q$一个闭环流程，接下来可以通过迭代法来做完成。针对上述例子，我们假设初始值为$\\pi^{0}=0.5，p^{0}=0.5，q^{0}=0.5$，因为对$y_i=1$和$y_i=0$均有$\\mu_j^{1}=0.5$，利用迭代公式计算得到$\\pi^{1}=0.5，p^{1}=0.6，q^{1}=0.6$，继续迭代得到最终的参数：\n",
    "$$\\widehat{\\pi^{0}}=0.5，\\widehat{p^{0}}=0.6，\\widehat{q^{0}}=0.6$$\n",
    "如果一开始初始值选择为：$\\pi^{0}=0.4，p^{0}=0.6，q^{0}=0.7$，那么得到的模型参数的极大似然估计是$\\widehat{\\pi}=0.4064，\\widehat{p}=0.5368，\\widehat{q}=0.6432$，这说明EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。\n",
    "\n",
    "这个例子中你只观察到了硬币抛完的结果，并不了解A硬币抛完之后，是选择了B硬币抛还是C硬币抛，这时候概率模型就存在着隐含变量！"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EM算法 \n",
    "输入：观测变量数据Y，隐变量数据Z，联合分布$P(Y,Z|\\theta)$，条件分布$P(Z|Y,\\theta)$；\n",
    "输出：模型参数$\\theta$\n",
    "- (1)选择参数的初值$\\theta^0$，开始迭代\n",
    "- (2) E步：记$\\theta^i$为第i次迭代参数$\\theta$的估计值，在第i+1次迭代的E步，计算\n",
    "$$\n",
    "\\begin{aligned}\n",
    "Q(\\theta,\\theta^i)&=E_{Z}[logP(Y,Z|\\theta)|Y,\\theta^i]\\\\\n",
    "&=\\sum_{Z}logP(Y,Z|\\theta)P(Z|Y,\\theta^i)\n",
    "\\end{aligned}\n",
    "$$\n",
    "这里，$P(Z|Y,\\theta^i)$是在给定观测数据Y和当前的参数估计$\\theta^i$下隐变量数据Z的条件概率分布；\n",
    "\n",
    "- (3) M步：求使$Q(\\theta,\\theta^i)$极大化的$\\theta$，确定第i+1次迭代的参数的估计值$\\theta^{i+1}$，\n",
    "$$\n",
    "\\theta^{i+1}=arg \\max \\limits_{\\theta}Q(\\theta,\\theta^{i})\n",
    "$$\n",
    "$Q(\\theta,\\theta^{i})$是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的。\n",
    "- (4) 重复第(2)步和第(3)步，直到收敛，收敛条件：\n",
    "$$\n",
    "|| \\theta^{i+1}-\\theta^{i} || < \\varepsilon_1\n",
    "$$\n",
    "或者：\n",
    "$$\n",
    "||Q(\\theta^{i+1},\\theta^{i})-Q(\\theta^{i},\\theta^{i})|| <\\varepsilon_2\n",
    "$$\n",
    "收敛迭代就结束了。我们来拆解一下这个M步骤，"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 推导逼近\n",
    "主要讲解Jensen不等式，这个公式在推导和收敛都用到，主要是如下的结论：\n",
    "- $f(x)$是凸函数\n",
    "$$\n",
    "f(E(X)) \\le E(f(x))\n",
    "$$\n",
    "- $f(x)$是凹函数\n",
    "$$\n",
    "f(E(X)) \\ge E(f(x))\n",
    "$$\n",
    "\n",
    "推导出Em算法可以近似实现对观测数据的极大似然估计的办法是找到E步骤的下界，让下届最大，通过逼近的方式实现对观测数据的最大似然估计。统计学习基础中采用的是相减方式，我们来看下具体的步骤。\n",
    "- 增加隐藏变量\n",
    "$$\n",
    "L(\\theta)=\\sum_{Z}logP(Y|Z,\\theta)P(Z,\\theta)\n",
    "$$\n",
    "则$L(\\theta)-L(\\theta^{i})$为：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "L(\\theta)-L(\\theta^{i})=log(\\sum_{Z} P(Y|Z,\\theta^i)\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Y|Z,\\theta^i)})-L(\\theta^{i})\\\\\n",
    "\\ge \\sum_{Z} P(Y|Z,\\theta^i)log(\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Y|Z,\\theta^i)})-L(\\theta^{i})\n",
    "\\end{aligned}\n",
    "$$\n",
    "$\\ge$这一个步骤就是采用了凹函数的Jensen不等式做转换。因为$Z$是隐藏变量，所以有$\\sum_{Z} P(Y|Z,\\theta^i)==1，P(Y|Z,\\theta^i)>0$，于是继续变：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "L(\\theta)-L(\\theta^{i})&=log(\\sum_{Z} P(Y|Z,\\theta^i)\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Y|Z,\\theta^i)})-L(\\theta^{i})\\\\\n",
    "&\\ge \\sum_{Z} P(Z|Y,\\theta^i)log(\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Z|Y,\\theta^i)})-L(\\theta^{i})\\\\\n",
    "&=\\sum_{Z} P(Z|Y,\\theta^i)log(\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Z|Y,\\theta^i)})-\\sum_{Z} P(Z|Y,\\theta^i)L(\\theta^{i})\\\\\n",
    "&= \\sum_{Z} P(Z|Y,\\theta^i)log(\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Z|Y,\\theta^i) (P(Y|\\theta^{i})}) \\\\\n",
    "& \\ge0\n",
    "\\end{aligned}\n",
    "$$\n",
    "也就是：$L(\\theta)\\ge L(\\theta^{i})+ \\sum_{Z} P(Z|Y,\\theta^i)log(\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Y|Z,\\theta^i) L(\\theta^{i})})$，有下界，最大化下界，来得到近似值。这里有一个细节：$P(Y|Z,\\theta^i)$ 变为$P(Z|Y,\\theta^i)$？如果要满足Jensen不等式的等号，则有：\n",
    "$$\n",
    "\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Y|Z,\\theta^i)} = c\n",
    "$$\n",
    "c为一个常数，而$\\sum_{Z}P(Y|Z,\\theta^i)=1$则：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\sum_{Z}P(Y|Z,\\theta)P(Z,\\theta)= c\\sum_{Z}P(Y|Z,\\theta^i)&=c\\\\\n",
    "&=\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Y|Z,\\theta^i)}\\\\\n",
    "P(Y|Z,\\theta)=\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{\\sum_{Z}P(Y|Z,\\theta)P(Z,\\theta)}=\\frac{P(Y,Z,\\theta)}{P(Y,\\theta)}=P(Z|Y,\\theta)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "大家是不是很奇怪$P(Y|Z,\\theta)P(Z,\\theta)$加上$\\sum$之后等于什么，其实有的博客这里使用$P(Z,\\theta) = P(Y^i,Z^i,\\theta^i)$来替代$P(Y|Z,\\theta)$参与计算，这样$\\sum_{Z}P(Y^i,Z^i,\\theta^i)$，这样就方便理解来了。\n",
    "\n",
    "于是最大化如下：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\theta^{i+1}&=arg \\max_{\\theta}\\sum_{Z} P(Z|Y,\\theta^i)log(\\frac{P(Y|Z,\\theta)P(Z,\\theta)}{P(Z|Y,\\theta^i)})\\\\\n",
    "&=arg \\max_{\\theta}\\sum_{Z} P(Z|Y,\\theta^i)log(P(Y|Z,\\theta)P(Z,\\theta))\\\\\n",
    "& =arg \\max_{\\theta}\\sum_{Z} P(Z|Y,\\theta^i)log(P(Y,Z|\\theta))\\\\ \n",
    "&=arg \\max_{\\theta}Q(\\theta,\\theta^i)\n",
    "\\end{aligned} \n",
    "$$\n",
    "其中$log$分母提出来是关于$Z$的$\\sum_{Z} P(Z|Y,\\theta^i)logP(Z|Y,\\theta^i)$，可以去掉。当然也有博客写的形式是：\n",
    "$$\n",
    "arg \\max_{\\theta}\\sum_{i=1}^{M}\\sum_{Z^{i}} P(Z^{i}|Y^{i},\\theta^i)log(P(Y^{i},Z^{i};\\theta))\\\\ \n",
    "$$\n",
    "形式其实一样，表示的不一样而已。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 证明收敛\n",
    "我们知道已知观测数据的似然函数是$P(Y,\\theta)$，对数似然函数为：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "L()=\\sum_{i=1}^{M}logP(y^{i},\\theta) &=\\sum_{i=1}^{M}log(\\frac{P(y^i,Z|\\theta)}{P(Z|y^i,\\theta)})\\\\\n",
    "&=\\sum_{i=1}^{M}logP(y^i,Z|\\theta) - \\sum_{i=1}^{M}logP(Z|y^i,\\theta)\n",
    "\\end{aligned} \n",
    "$$\n",
    "要证明收敛，就证明单调递增，$\\sum_{i=1}^{M}logP(y^{i},\\theta^{j+1})>\\sum_{i=1}^{M}logP(y^{i},\\theta^{j})$\n",
    "由上文知道：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "Q(\\theta,\\theta^i)&=\\sum_{Z}logP(Y,Z|\\theta)P(Z|Y,\\theta^i)\\\\\n",
    "&=\\sum_{i=1}^{M}\\sum_{Z^j}logP(y^i,Z^j|\\theta)P(Z^j|y^i,\\theta^i)\n",
    "\\end{aligned}\n",
    "$$\n",
    "我们构造一个函数$H$，让他等于：\n",
    "$$\n",
    "H(\\theta,\\theta^{i})=\\sum_{i=1}^{M}\\sum_{Z^j}log(P(Z|y^i,\\theta)P(Z|y^i,\\theta^i))\n",
    "$$\n",
    "让$Q(\\theta,\\theta^i)-H(\\theta,\\theta^{i})$：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "Q(\\theta,\\theta^i)-H(\\theta,\\theta^{i})&=\\sum_{i=1}^{M}\\sum_{Z^j}logP(y^i,Z^j|\\theta)P(Z^j|y^i,\\theta^i) - \\sum_{i=1}^{M}\\sum_{Z^j}log(P(Z^j|y^i,\\theta)P(Z^j|y^i,\\theta^i)) \\\\\n",
    "&=\\sum_{i=1}^{M}\\sum_{Z^j}log\\bigg(P(y^i,Z^j|\\theta)-P(Z^j|y^i,\\theta)\\bigg) \\\\\n",
    "&=\\sum_{i=1}^{M}logP(y^{i},\\theta) \n",
    "\\end{aligned}\n",
    "$$所以：\n",
    "$$\n",
    "\\sum_{i=1}^{M}logP(y^{i},\\theta^{j+1})-\\sum_{i=1}^{M}logP(y^{i},\\theta^{j}) \\\\\n",
    "= Q(\\theta^{i+1},\\theta^i)-H(\\theta^{i+1},\\theta^{i}) - (Q(\\theta^{i},\\theta^{i})-H(\\theta^{i},\\theta^{i}))\\\\\n",
    "= Q(\\theta^{i+1},\\theta^i)- Q(\\theta^{i},\\theta^{i}) -( H(\\theta^{i+1},\\theta^{i}) - H(\\theta^{i},\\theta^{i}))\n",
    "$$\n",
    "该公式左边已经被证明是大于0，证明右边：$H(\\theta^{i+1},\\theta^{i}) - H(\\theta^{i},\\theta^{i})<0$：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "H(\\theta^{i+1},\\theta^{i}) - H(\\theta^{i},\\theta^{i}) &=\\sum_{Z^j}\\bigg(log(\\frac{P(Z^j|Y,\\theta^{i+1})}{P(Z^j|Y,\\theta^i)}) \\bigg)P(Z^j|Y,\\theta^i) \\\\\n",
    "&=log\\bigg(\\sum_{Z^j}\\frac{P(Z^j|Y,\\theta^{i+1})}{P(Z^j|Y,\\theta^i)}P(Z^j|Y,\\theta^i) \\bigg)\\\\\n",
    "&=logP(Z|Y,\\theta^{i+1})=log1=0\n",
    "\\end{aligned}\n",
    "$$\n",
    "其中不等式是由于Jensen不等式，由此证明了$\\sum_{i=1}^{M}logP(y^{i},\\theta^{j+1})>\\sum_{i=1}^{M}logP(y^{i},\\theta^{j})$，证明了EM算法的收敛性。但不能保证是全局最优，只能保证局部最优。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 高斯混合分布\n",
    "EM算法的一个重要应用场景就是高斯混合模型的参数估计。高斯混合模型就是由多个高斯模型组合在一起的混合模型（可以理解为多个高斯分布函数的线性组合，理论上高斯混合模型是可以拟合任意类型的分布），例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的：\n",
    "![](https://imgconvert.csdnimg.cn/aHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL0tsYXVzemhhby9waWN0dXJlL21hc3Rlci9waWN0dXJlL2NvbW1vbi8lRTklQUIlOTglRTYlOTYlQUYlRTYlQjclQjclRTUlOTAlODglRTUlOEQlOTUlRTUlODglODYlRTUlQjglODMucG5n?x-oss-process=image/format,png)\n",
    "\n",
    "两个高斯模型可以拟合数据集，如图所示：\n",
    "![](https://imgconvert.csdnimg.cn/aHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL0tsYXVzemhhby9waWN0dXJlL21hc3Rlci9waWN0dXJlL2NvbW1vbi8lRTklQUIlOTglRTYlOTYlQUYlRTYlQjclQjclRTUlOTAlODglRTUlODglODYlRTUlQjglODMtJUU1JThGJThDJUU1JTg4JTg2JUU1JUI4JTgzLnBuZw?x-oss-process=image/format,png)\n",
    "\n",
    "如果有多个高斯模型，公式表示为：\n",
    "$$\n",
    "P(y|\\theta)=\\sum_{k=1}^{K}a_k\\phi(y|\\theta_{k}) \\\\\n",
    "\\phi(y|\\theta_{k})=\\frac{1}{\\sqrt{2\\pi}\\delta_{k}}exp(-\\frac{(y-\\mu_{k})^2}{2 \\delta_{k}^{2}}) \\\\\n",
    "a_k>0,\\sum a_k =1\n",
    "$$\n",
    "$\\phi(y|\\theta_{k})$表示为第k个高斯分布密度模型，定义如上，其中$a_k$表示被选中的概率。在本次模型$P(y|\\theta)$中，观测数据是已知的，而观测数据具体来自哪个模型是未知的，有点像之前提过的三硬币模型，我们来对比一下，A硬币就像是概率$a_k$，用来表明具体的模型，而B、C硬币就是具体的模型，只不过这里有很多个模型，不仅仅是B、C这两个模型。我们用$\\gamma_{jk}$来表示，则：\n",
    "$$\n",
    "\\gamma_{jk} =\n",
    "\\begin{cases}\n",
    "1& \\text{第j个观测数据来源于第k个模型}\\\\\n",
    "0& \\text{否则}\n",
    "\\end{cases}\n",
    "$$\n",
    "所以一个观测数据$y_j$的隐藏数据$(\\gamma_{j1},\\gamma_{j2},...,\\gamma_{jk})$，那么完全似然函数就是：\n",
    "\n",
    "$$\n",
    "P(y,\\gamma|\\theta)= \\prod_{k=1}^{K}\\prod_{j=1}^{N}[a_{k}\\phi(y|\\theta_{k})]^{\\gamma_{jk}}\n",
    "$$\n",
    "\n",
    "取对数之后等于：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "log(P(y,\\gamma|\\theta))&=log( \\prod_{k=1}^{K}\\prod_{j=1}^{N}[a_{k}\\phi(y|\\theta_{k})]^{\\gamma_{jk}})\\\\\n",
    "&=\\sum_{K}^{k=1}\\bigg(\\sum_{j=1}^{N}(\\gamma_{jk}) log(a_k)+\\sum_{j=1}^{N}( \\gamma_{jk})\\bigg[log(\\frac{1}{\\sqrt{2\\pi}})-log(\\delta_{k})-\\frac{(y_i-\\mu_{k})^2}{2 \\delta_{k}^{2}}\\bigg]\\bigg)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "- E 步 ：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "Q(\\theta.\\theta^i) &= E[log(P(y,\\gamma|\\theta))]\\\\\n",
    "&=\\sum_{K}^{k=1}\\bigg(\\sum_{j=1}^{N}(E\\gamma_{jk}) log(a_k)+\\sum_{j=1}^{N}(E\\gamma_{jk})\\bigg[log(\\frac{1}{\\sqrt{2\\pi}})-log(\\delta_{k})-\\frac{(y_i-\\mu_{k})^2}{2 \\delta_{k}^{2}}\\bigg]\\bigg)\n",
    "\\end{aligned}\n",
    "$$\n",
    "其中我们定义$\\hat{\\gamma_{jk}}$：\n",
    "$$\n",
    "\\hat{\\gamma_{jk}} = E(\\gamma_{jk}|y,\\theta)=\\frac{a_k\\phi(y_i|\\theta_{k})}{\\sum_{k=1}^{K}a_k\\phi(y_i|\\theta_{k}) }\\\\\n",
    "j=1,2,..,N；k=1,2,...,K\\\\\n",
    "n_k=\\sum_{j=i}^{N}E\\gamma_{jk}\n",
    "$$\n",
    "于是化简得到：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "Q(\\theta.\\theta^i) &= \\sum_{K}^{k=1}\\bigg(n_k log(a_k)+\\sum_{j=1}^{N}(E\\gamma_{jk})\\bigg[log(\\frac{1}{\\sqrt{2\\pi}})-log(\\delta_{k})-\\frac{(y_i-\\mu_{k})^2}{2 \\delta_{k}^{2}}\\bigg]\\bigg)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "E 步 在代码设计上只有$\\hat{\\gamma_{jk}}$有用，用于M步的计算。\n",
    "\n",
    "\n",
    "- M步，\n",
    "$$\n",
    "\\theta^{i+1}=arg \\max_{\\theta}Q(\\theta,\\theta^i)\n",
    "$$\n",
    "对$Q(\\theta,\\theta^i)$求导，得到每个未知量的偏导，使其偏导等于0，求解得到：\n",
    "$$\n",
    "\\hat{\\mu_k}=\\frac{\\sum_{j=1}^{N}\\hat{\\gamma_{jk}}y_i}{\\sum_{j=1}^{N}\\hat{\\gamma_{jk}}}\n",
    "\\\\\n",
    "\\\\\n",
    "\\hat{\\delta_k}=\\frac{\\sum_{j=1}^{N}\\hat{\\gamma_{jk}}(y_i-\\mu_k)^2}{\\sum_{j=1}^{N}\\hat{\\gamma_{jk}}}\n",
    "\\\\\n",
    "\\\\\n",
    "\\\\\n",
    "\\hat{a_k}=\\frac{\\sum_{j=1}^{N}\\hat{\\gamma_{jk}} }{N}\n",
    "$$\n",
    "给一个初始值，来回迭代就可以求得值内容。这一块主要用到了$Q(\\theta.\\theta^i)$的导数，并且用到了E步的$\\hat{\\gamma_{jk}}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 混合高斯分布模型\n",
    "EM算法更多是一种思想，用概率来解决问题的一种方法，具体的代码看自己选用模型，所以并没有通用的模型，本此代码主要是讲解混合高斯分布模型的\n",
    "\n",
    "这其中的M步 完全按照了 公式来计算。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "import math\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "数据集：伪造数据集（两个高斯分布混合）\n",
    "数据集长度：1000\n",
    "------------------------------\n",
    "运行结果：\n",
    "----------------------------\n",
    "the Parameters set is:\n",
    "alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0\n",
    "----------------------------\n",
    "the Parameters predict is:\n",
    "alpha0:0.4, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9\n",
    "----------------------------\n",
    "'''\n",
    "\n",
    "def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):\n",
    "    '''\n",
    "    初始化数据集\n",
    "    这里通过服从高斯分布的随机函数来伪造数据集\n",
    "    :param mu0: 高斯0的均值\n",
    "    :param sigma0: 高斯0的方差\n",
    "    :param mu1: 高斯1的均值\n",
    "    :param sigma1: 高斯1的方差\n",
    "    :param alpha0: 高斯0的系数\n",
    "    :param alpha1: 高斯1的系数\n",
    "    :return: 混合了两个高斯分布的数据\n",
    "    '''\n",
    "    # 定义数据集长度为1000\n",
    "    length = 1000\n",
    "\n",
    "    # 初始化第一个高斯分布，生成数据，数据长度为length * alpha系数，以此来\n",
    "    # 满足alpha的作用\n",
    "    data0 = np.random.normal(mu0, sigma0, int(length * alpha0))\n",
    "    # 第二个高斯分布的数据\n",
    "    data1 = np.random.normal(mu1, sigma1, int(length * alpha1))\n",
    "\n",
    "    # 初始化总数据集\n",
    "    # 两个高斯分布的数据混合后会放在该数据集中返回\n",
    "    dataSet = []\n",
    "    # 将第一个数据集的内容添加进去\n",
    "    dataSet.extend(data0)\n",
    "    # 添加第二个数据集的数据\n",
    "    dataSet.extend(data1)\n",
    "    # 对总的数据集进行打乱（其实不打乱也没事，只不过打乱一下直观上让人感觉已经混合了\n",
    "    # 读者可以将下面这句话屏蔽以后看看效果是否有差别）\n",
    "    random.shuffle(dataSet)\n",
    "\n",
    "    #返回伪造好的数据集\n",
    "    return dataSet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 高斯分布公式，没有什么特殊的\n",
    "def calcGauss(dataSetArr, mu, sigmod):\n",
    "    '''\n",
    "    根据高斯密度函数计算值\n",
    "    依据：“9.3.1 高斯混合模型” 式9.25\n",
    "    注：在公式中y是一个实数，但是在EM算法中(见算法9.2的E步)，需要对每个j\n",
    "    都求一次yjk，在本实例中有1000个可观测数据，因此需要计算1000次。考虑到\n",
    "    在E步时进行1000次高斯计算，程序上比较不简洁，因此这里的y是向量，在numpy\n",
    "    的exp中如果exp内部值为向量，则对向量中每个值进行exp，输出仍是向量的形式。\n",
    "    所以使用向量的形式1次计算即可将所有计算结果得出，程序上较为简洁\n",
    "    \n",
    "    :param dataSetArr: 可观测数据集\n",
    "    :param mu: 均值\n",
    "    :param sigmod: 方差\n",
    "    :return: 整个可观测数据集的高斯分布密度（向量形式）\n",
    "    '''\n",
    "    # 计算过程就是依据式9.25写的，没有别的花样\n",
    "    result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))\n",
    "    # 返回结果\n",
    "    return result\n",
    "\n",
    "\n",
    "def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):\n",
    "    '''\n",
    "    EM算法中的E步\n",
    "    依据当前模型参数，计算分模型k对观数据y的响应度\n",
    "    :param dataSetArr: 可观测数据y\n",
    "    :param alpha0: 高斯模型0的系数\n",
    "    :param mu0: 高斯模型0的均值\n",
    "    :param sigmod0: 高斯模型0的方差\n",
    "    :param alpha1: 高斯模型1的系数\n",
    "    :param mu1: 高斯模型1的均值\n",
    "    :param sigmod1: 高斯模型1的方差\n",
    "    :return: 两个模型各自的响应度\n",
    "    '''\n",
    "    # 计算y0的响应度\n",
    "    # 先计算模型0的响应度的分子\n",
    "    gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)\n",
    "    #print(\"gamma0=\",gamma0.shape) # 1000, 维向量\n",
    "    # 模型1响应度的分子\n",
    "    gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)\n",
    "\n",
    "    # 两者相加为E步中的分布\n",
    "    sum = gamma0 + gamma1\n",
    "    # 各自相除，得到两个模型的响应度\n",
    "    gamma0 = gamma0 / sum\n",
    "    gamma1 = gamma1 / sum\n",
    "\n",
    "    # 返回两个模型响应度\n",
    "    return gamma0, gamma1\n",
    "\n",
    "def M_step(muo, mu1, gamma0, gamma1, dataSetArr):\n",
    "    # 依据算法9.2计算各个值\n",
    "    # 这里没什么花样，对照书本公式看看这里就好了\n",
    "    \n",
    "    # np.dot 点积：[1,2] [2,3] = [2,6]\n",
    "    mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)\n",
    "    mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)\n",
    "\n",
    "    # math.sqrt  平方根 \n",
    "    sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0))\n",
    "    sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1))\n",
    "\n",
    "    alpha0_new = np.sum(gamma0) / len(gamma0)\n",
    "    alpha1_new = np.sum(gamma1) / len(gamma1)\n",
    "\n",
    "    # 将更新的值返回\n",
    "    return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new\n",
    "\n",
    "\n",
    "## 训练主函数\n",
    "def EM_Train(dataSetList, iter=500):\n",
    "    '''\n",
    "    根据EM算法进行参数估计\n",
    "    算法依据“9.3.2 高斯混合模型参数估计的EM算法” 算法9.2\n",
    "    :param dataSetList:数据集（可观测数据）\n",
    "    :param iter: 迭代次数\n",
    "    :return: 估计的参数\n",
    "    '''\n",
    "    # 将可观测数据y转换为数组形式，主要是为了方便后续运算\n",
    "    dataSetArr = np.array(dataSetList)\n",
    "\n",
    "    # 步骤1：对参数取初值，开始迭代\n",
    "    alpha0 = 0.5\n",
    "    mu0 = 0\n",
    "    sigmod0 = 1\n",
    "    alpha1 = 0.5\n",
    "    mu1 = 1\n",
    "    sigmod1 = 1\n",
    "\n",
    "    # 开始迭代\n",
    "    step = 0\n",
    "    while (step < iter):\n",
    "        # 每次进入一次迭代后迭代次数加1\n",
    "        step += 1\n",
    "        # 步骤2：E步：依据当前模型参数，计算分模型k对观测数据y的响应度\n",
    "        gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)\n",
    "        # 步骤3：M步\n",
    "        mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)\n",
    "\n",
    "    # 迭代结束后将更新后的各参数返回\n",
    "    return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------\n",
      "the Parameters set is:\n",
      "alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0\n",
      "----------------------------\n",
      "the Parameters predict is:\n",
      "alpha0:0.4, mu0:0.6, sigmod0:-1.8, alpha1:0.8, mu1:0.7, sigmod1:0.9\n",
      "----------------------------\n",
      "time span: 0.26200032234191895\n"
     ]
    }
   ],
   "source": [
    "if __name__ == '__main__':\n",
    "    start = time.time()\n",
    "\n",
    "    # 设置两个高斯模型进行混合，这里是初始化两个模型各自的参数\n",
    "    # 见“9.3 EM算法在高斯混合模型学习中的应用”\n",
    "    # alpha是“9.3.1 高斯混合模型” 定义9.2中的系数α\n",
    "    # mu0是均值μ\n",
    "    # sigmod是方差σ\n",
    "    # 在设置上两个alpha的和必须为1，其他没有什么具体要求，符合高斯定义就可以\n",
    "    \n",
    "    alpha0 = 0.3  # 系数α\n",
    "    mu0 = -2  # 均值μ\n",
    "    sigmod0 = 0.5  # 方差σ\n",
    "\n",
    "    alpha1 = 0.7  # 系数α\n",
    "    mu1 = 0.5  # 均值μ\n",
    "    sigmod1 = 1  # 方差σ\n",
    "\n",
    "    # 初始化数据集\n",
    "    dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)\n",
    "\n",
    "    #打印设置的参数\n",
    "    print('---------------------------')\n",
    "    print('the Parameters set is:')\n",
    "    print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (\n",
    "        alpha0, alpha1, mu0, mu1, sigmod0, sigmod1\n",
    "    ))\n",
    "\n",
    "    # 开始EM算法，进行参数估计\n",
    "    alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)\n",
    "\n",
    "    # 打印参数预测结果\n",
    "    print('----------------------------')\n",
    "    print('the Parameters predict is:')\n",
    "    print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (\n",
    "        alpha0, alpha1, mu0, mu1, sigmod0, sigmod1\n",
    "    ))\n",
    "\n",
    "    # 打印时间\n",
    "    print('----------------------------')\n",
    "    print('time span:', time.time() - start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## E步主要计算内容\n",
    "其中我们定义$\\hat{\\gamma_{jk}}$：\n",
    "$$\n",
    "\\hat{\\gamma_{jk}} = E(\\gamma_{jk}|y,\\theta)=\\frac{a_k\\phi(y_i|\\theta_{k})}{\\sum_{k=1}^{K}a_k\\phi(y_i|\\theta_{k}) }\\\\\n",
    "j=1,2,..,N；k=1,2,...,K\\\\\n",
    "n_k=\\sum_{j=i}^{N}E\\gamma_{jk}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## M步 主要计算内容\n",
    "这一步骤主要是对Q函数求导后的数据进行计算，利用了 E 步 的$\\hat{\\gamma_{jk}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import copy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    " \n",
    "#生成随机数据，4个高斯模型\n",
    "def generate_data(sigma,N,mu1,mu2,mu3,mu4,alpha):\n",
    "    global X                  #可观测数据集\n",
    "    X = np.zeros((N, 2))       # 初始化X，2行N列。2维数据，N个样本\n",
    "    X=np.matrix(X)\n",
    "    global mu                 #随机初始化mu1，mu2，mu3，mu4\n",
    "    mu = np.random.random((4,2))\n",
    "    mu=np.matrix(mu)\n",
    "    global excep              #期望第i个样本属于第j个模型的概率的期望\n",
    "    excep=np.zeros((N,4))\n",
    "    global alpha_             #初始化混合项系数\n",
    "    alpha_=[0.25,0.25,0.25,0.25]\n",
    "    for i in range(N):\n",
    "        if np.random.random(1) < 0.1:  # 生成0-1之间随机数\n",
    "            X[i,:]  = np.random.multivariate_normal(mu1, sigma, 1)     #用第一个高斯模型生成2维数据\n",
    "        elif 0.1 <= np.random.random(1) < 0.3:\n",
    "            X[i,:] = np.random.multivariate_normal(mu2, sigma, 1)      #用第二个高斯模型生成2维数据\n",
    "        elif 0.3 <= np.random.random(1) < 0.6:\n",
    "            X[i,:] = np.random.multivariate_normal(mu3, sigma, 1)      #用第三个高斯模型生成2维数据\n",
    "        else:\n",
    "            X[i,:] = np.random.multivariate_normal(mu4, sigma, 1)      #用第四个高斯模型生成2维数据\n",
    " \n",
    "    print(\"可观测数据：\\n\",X)       #输出可观测样本\n",
    "    print(\"初始化的mu1，mu2，mu3，mu4：\",mu)      #输出初始化的mu\n",
    "\n",
    "\n",
    "# E 期望\n",
    "#  \\hat{\\gamma_{jk}}\n",
    "def e_step(sigma,k,N):\n",
    "    global X\n",
    "    global mu\n",
    "    global excep\n",
    "    global alpha_\n",
    "    for i in range(N):\n",
    "        denom=0\n",
    "        for j in range(0,k):\n",
    "            #  sigma.I 表示矩阵的逆矩阵\n",
    "            # np.transpose ：矩阵转置   np.linalg.det():矩阵求行列式\n",
    "            denom += alpha_[j]*  math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))  /np.sqrt(np.linalg.det(sigma))       #分母\n",
    "        for j in range(0,k):\n",
    "            numer = math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/np.sqrt(np.linalg.det(sigma))        #分子\n",
    "            excep[i,j]=alpha_[j]*numer/denom      #求期望\n",
    "    print(\"隐藏变量：\\n\",excep)\n",
    "\n",
    "    \n",
    "def m_step(k,N):\n",
    "    global excep\n",
    "    global X\n",
    "    global alpha_\n",
    "    for j in range(0,k):\n",
    "        denom=0   #分母\n",
    "        numer=0   #分子\n",
    "        for i in range(N):\n",
    "            numer += excep[i,j]*X[i,:]\n",
    "            denom += excep[i,j]\n",
    "        mu[j,:] = numer/denom    #求均值\n",
    "        alpha_[j]=denom/N        #求混合项系数\n",
    "\n",
    "        #     #可视化结果\n",
    "def plotShow():\n",
    "    # 画生成的原始数据\n",
    "    plt.subplot(221)\n",
    "    plt.scatter(X[:,0].tolist(), X[:,1].tolist(),c='b',s=25,alpha=0.4,marker='o')    #T散点颜色，s散点大小，alpha透明度，marker散点形状\n",
    "    plt.title('random generated data')\n",
    "    #画分类好的数据\n",
    "    plt.subplot(222)\n",
    "    plt.title('classified data through EM')\n",
    "    order=np.zeros(N)\n",
    "    color=['b','r','k','y']\n",
    "    for i in range(N):\n",
    "        for j in range(k):\n",
    "            if excep[i,j]==max(excep[i,:]):\n",
    "                order[i]=j     #选出X[i,:]属于第几个高斯模型\n",
    "            probility[i] += alpha_[int(order[i])]*math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/(np.sqrt(np.linalg.det(sigma))*2*np.pi)    #计算混合高斯分布\n",
    "        plt.scatter(X[i, 0], X[i, 1], c=color[int(order[i])], s=25, alpha=0.4, marker='o')      #绘制分类后的散点图\n",
    "    #绘制三维图像\n",
    "    ax = plt.subplot(223, projection='3d')\n",
    "    plt.title('3d view')\n",
    "    for i in range(N):\n",
    "        ax.scatter(X[i, 0], X[i, 1], probility[i], c=color[int(order[i])])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "可观测数据：\n",
      " [[27.83586245 38.22407571]\n",
      " [14.36625986 19.67437442]\n",
      " [24.60609231 23.73534404]\n",
      " [16.7951181  36.01232204]\n",
      " [26.47662418 16.70777218]\n",
      " [51.47497314 13.10622654]\n",
      " [39.317069   20.62520832]\n",
      " [19.78348463 28.92277082]\n",
      " [27.87073258 17.17796032]\n",
      " [29.78511733 38.75742229]\n",
      " [51.41115437 21.00448716]\n",
      " [36.87306609 24.48003687]\n",
      " [16.31532846 16.51489081]\n",
      " [45.21454785 13.69380145]\n",
      " [40.42445226 14.06947943]\n",
      " [44.5990409  20.4189464 ]\n",
      " [42.482583    1.99874287]\n",
      " [30.94584559 36.90250786]\n",
      " [ 9.40084917 28.26212545]\n",
      " [46.40015649 13.74067493]\n",
      " [27.05400361 39.58684097]\n",
      " [44.78350628 17.65849471]\n",
      " [40.13141565 14.42945166]\n",
      " [13.29592491 21.6924786 ]\n",
      " [56.72572468  9.37636317]\n",
      " [42.79571057 19.59989619]\n",
      " [44.43058709  8.11809346]\n",
      " [ 8.41883723 22.42134984]\n",
      " [ 2.30194905 34.86979563]\n",
      " [29.02601979 30.83736503]\n",
      " [42.34076607 20.44169753]\n",
      " [ 9.12726519 31.58709971]\n",
      " [34.43623505 17.19609801]\n",
      " [50.69173891 11.37430969]\n",
      " [47.05370416 17.03456568]\n",
      " [39.38806226 20.19467327]\n",
      " [16.94771874 16.67321045]\n",
      " [19.80650541 24.22988707]\n",
      " [25.50054716 33.51271549]\n",
      " [ 7.83222534 30.90246427]\n",
      " [45.45475008  8.09911822]\n",
      " [ 3.32716342 35.7496015 ]\n",
      " [37.20659713 13.62374041]\n",
      " [30.4812938  40.70390448]\n",
      " [40.37749213 12.62125773]\n",
      " [41.49907331 13.33930952]\n",
      " [ 9.00012529 31.75836259]\n",
      " [46.0189854  18.56505336]\n",
      " [50.14228211  7.96583182]\n",
      " [51.43954113 20.25449608]\n",
      " [20.51458706 32.3994836 ]\n",
      " [41.71503817  7.61408459]\n",
      " [45.20623519 37.86509752]\n",
      " [28.83490553 40.69030451]\n",
      " [28.84252313 38.55073538]\n",
      " [50.49391528 18.59097389]\n",
      " [18.82244005 20.5566116 ]\n",
      " [45.2483832   9.69162326]\n",
      " [32.68733826  9.42413435]\n",
      " [13.09742812 18.26478552]\n",
      " [39.34671103 18.50622696]\n",
      " [30.08543924 33.24647564]\n",
      " [49.85453826  3.00210314]\n",
      " [12.92815575 33.60088853]\n",
      " [24.01736904 23.45685367]\n",
      " [14.20669024 11.38883184]\n",
      " [48.99216386 15.78809094]\n",
      " [50.80161787  6.81695344]\n",
      " [ 6.17129427 35.32445819]\n",
      " [49.61416008 20.3279607 ]\n",
      " [41.20556142 24.08539969]\n",
      " [36.53178968 46.84189442]\n",
      " [52.08160516 15.25012141]\n",
      " [-0.85558581 30.69761062]\n",
      " [50.67042641 18.05017014]\n",
      " [49.23380909 24.75131372]\n",
      " [43.89387202  3.64089196]\n",
      " [52.96119982  9.22461525]\n",
      " [56.1972556  12.74935673]\n",
      " [46.43947201 19.06824739]\n",
      " [ 9.85683016 35.88448577]\n",
      " [24.21872348 20.37951855]\n",
      " [11.66161124 36.2232858 ]\n",
      " [46.2240621  18.46000115]\n",
      " [ 8.77376786 35.68975593]\n",
      " [ 4.03923919 36.47108537]\n",
      " [42.31416487 24.48839236]\n",
      " [43.57418002 11.03194241]\n",
      " [25.54031836 20.36223642]\n",
      " [43.37291411 19.60807727]\n",
      " [20.3781315  19.89139521]\n",
      " [40.45929215 16.20410483]\n",
      " [48.70914629 13.77714863]\n",
      " [24.27327304 38.65121008]\n",
      " [40.08677883 11.58070999]\n",
      " [24.10964775 36.84763744]\n",
      " [52.65783202 14.14041665]\n",
      " [ 1.44977753 30.62066858]\n",
      " [43.75798746 15.59303747]\n",
      " [15.37222662 18.45588884]\n",
      " [46.47258013 13.88919317]\n",
      " [44.22232092 12.95212251]\n",
      " [27.11998226 42.38657982]\n",
      " [41.90774788 16.86166716]\n",
      " [35.25947391 17.22575757]\n",
      " [ 0.84190096 32.63486972]\n",
      " [42.65819744 11.74118108]\n",
      " [ 8.26962447 33.63846321]\n",
      " [27.3772952  15.01254383]\n",
      " [37.58555751 12.04724541]\n",
      " [48.6360956  15.08563641]\n",
      " [ 1.22280721 36.5176344 ]\n",
      " [-2.03516614 33.33540002]\n",
      " [51.74347947 17.03723179]\n",
      " [11.67712313 16.83360469]\n",
      " [14.43369034 17.10468754]\n",
      " [44.99024895 13.83050679]\n",
      " [19.91269458 19.94589626]\n",
      " [32.06111971 29.90396701]\n",
      " [16.89942624 32.5743425 ]\n",
      " [36.79128436 20.059835  ]\n",
      " [42.13023181 24.37628091]\n",
      " [44.76459255  8.24460501]\n",
      " [18.10304029 18.14893112]\n",
      " [55.43784583  2.64425731]\n",
      " [45.77251166 18.75664651]\n",
      " [17.55660325 46.38741297]\n",
      " [ 9.84412627 37.01632286]\n",
      " [42.06530239 14.96542249]\n",
      " [-0.44792299 39.53764637]\n",
      " [43.59346593 11.55049217]\n",
      " [48.33827133 10.43941702]\n",
      " [46.69314219 16.13284812]\n",
      " [37.18683239 34.31605695]\n",
      " [37.69399903 19.86496582]\n",
      " [43.18795261  8.54896507]\n",
      " [45.25007406 23.26879671]\n",
      " [18.39254285 16.4594524 ]\n",
      " [54.45551139  8.08314286]\n",
      " [48.85805931 12.7935342 ]\n",
      " [28.36487701 45.80999821]\n",
      " [40.32439416 12.21443996]\n",
      " [ 5.54102391 36.44802424]\n",
      " [11.95152055 36.1625337 ]\n",
      " [37.51299858 12.00305227]\n",
      " [ 5.92617161 38.71806325]\n",
      " [23.28532635 34.26978638]\n",
      " [35.96348362 14.99084425]\n",
      " [41.87406558 16.43474897]\n",
      " [ 3.21217661 36.70757833]\n",
      " [27.36414372 40.36432639]\n",
      " [29.70448178 43.65291529]\n",
      " [31.40864825 19.00060567]\n",
      " [43.53437843 13.14585141]\n",
      " [19.50050805 19.78641395]\n",
      " [23.75768366 39.37069199]\n",
      " [40.85390311 15.18350218]\n",
      " [55.97603798 20.38022024]\n",
      " [41.20459654 22.44285663]\n",
      " [15.72483143 20.18865777]\n",
      " [ 7.09976542 25.25318266]\n",
      " [23.73458271 22.60728732]\n",
      " [44.65269778 21.9308162 ]\n",
      " [26.83763828 41.81800881]\n",
      " [47.44141999 14.20723093]\n",
      " [45.38314173 23.36699268]\n",
      " [46.30567427 26.63499278]\n",
      " [20.38915828 20.14797422]\n",
      " [20.07137706 12.28884043]\n",
      " [23.60257429 27.75115367]\n",
      " [29.06126738 32.53743365]\n",
      " [28.63281996 38.07111116]\n",
      " [48.63224748 12.35495171]\n",
      " [47.73065204 21.17417998]\n",
      " [41.94823644 13.96412374]\n",
      " [40.57364993  6.9461019 ]\n",
      " [30.58887075 38.54868104]\n",
      " [50.07953112 16.28051908]\n",
      " [34.16180569  6.60255072]\n",
      " [46.4092203  15.43645859]\n",
      " [12.90588366 24.20425889]\n",
      " [45.04005196 15.92640369]\n",
      " [-1.90378658 42.57711893]\n",
      " [30.51867385 35.65978162]\n",
      " [ 7.06937593 41.30247973]\n",
      " [16.21848423 22.57823424]\n",
      " [34.70469551 41.63859498]\n",
      " [-3.44636959 36.15251075]\n",
      " [47.06269576 22.40145536]\n",
      " [41.2063224   9.56669724]\n",
      " [ 6.62405736 32.25892335]\n",
      " [26.1016044  14.44488077]\n",
      " [44.81252704 14.42975033]\n",
      " [32.89356174 24.11344274]\n",
      " [32.7389753  31.857489  ]\n",
      " [42.46720618 16.97511204]\n",
      " [43.00284833 19.2662473 ]\n",
      " [21.16657366 23.05510421]\n",
      " [23.23976684 22.40216794]\n",
      " [29.22778928 46.80157791]\n",
      " [31.5771621  45.9657234 ]\n",
      " [16.3260274  42.80400325]\n",
      " [51.30782508 18.46277003]\n",
      " [15.00860037 25.58085553]\n",
      " [45.58438033 20.05478147]\n",
      " [36.56599006 17.04530652]\n",
      " [47.75784446  6.01278882]\n",
      " [40.70247319 18.18342598]\n",
      " [38.43220462  5.78224319]\n",
      " [50.49075371  7.49338251]\n",
      " [19.4047027  15.46216632]\n",
      " [42.93636988 19.34017643]\n",
      " [38.09858739 44.69362101]\n",
      " [42.47489889 11.58371963]\n",
      " [40.2218064  17.05735861]\n",
      " [54.84120491 16.95657372]\n",
      " [ 4.01703722 43.0009696 ]\n",
      " [ 7.7809159  38.83379979]\n",
      " [20.17178343 39.71047579]\n",
      " [17.34885016  7.67991286]\n",
      " [ 8.43707936 47.0362502 ]\n",
      " [43.84484023 20.47474344]\n",
      " [43.44093361 18.0195397 ]\n",
      " [39.38499471 10.35936637]\n",
      " [17.88646997 26.71876597]\n",
      " [47.36108528  8.34241341]\n",
      " [26.60302579 45.42210517]\n",
      " [ 5.08016004 42.54352087]\n",
      " [38.77839843 14.1463928 ]\n",
      " [46.56074267  9.33185451]\n",
      " [33.44036767 45.82645589]\n",
      " [36.31466859 19.07462415]\n",
      " [53.64527093 20.78719799]\n",
      " [51.36683289 18.27468071]\n",
      " [17.52537928 18.72744058]\n",
      " [23.41235699 22.53434668]\n",
      " [49.65293069 20.90399016]\n",
      " [34.74897857 19.89626669]\n",
      " [47.26753144 18.90977116]\n",
      " [43.04394107 25.22476244]\n",
      " [19.28043219 11.15785215]\n",
      " [45.42791407 19.01392978]\n",
      " [44.65428383 10.96126169]\n",
      " [44.83345978  6.8572207 ]\n",
      " [41.06044993 12.5306202 ]\n",
      " [36.87384999 19.92418527]\n",
      " [36.67782754 28.74630119]\n",
      " [-5.54671001 31.79746625]\n",
      " [27.78284626 36.53536253]\n",
      " [39.51285929 20.15624792]\n",
      " [15.24543924 27.59511125]\n",
      " [18.37860699 22.82786647]\n",
      " [49.01316178 19.01378331]\n",
      " [23.15429177 22.19790474]\n",
      " [19.71793047 23.41236615]\n",
      " [ 0.53042046 32.27864961]\n",
      " [19.8628091  13.85589626]\n",
      " [46.58355189 12.18125252]\n",
      " [46.69883918 18.04805034]\n",
      " [18.06668185 13.85382505]\n",
      " [43.93321513 12.75624356]\n",
      " [43.98688316  7.31350921]\n",
      " [38.70118051 34.41610495]\n",
      " [44.69839411 11.16574398]\n",
      " [43.77639576  8.12197135]\n",
      " [16.37947195 14.87600651]\n",
      " [38.6844825  41.3398049 ]\n",
      " [47.8082616  10.76492757]\n",
      " [ 8.68189161 50.62963987]\n",
      " [48.37878168 22.96682276]\n",
      " [37.995975   22.82467938]\n",
      " [47.44382325 11.44482419]\n",
      " [48.61339605 11.31256723]\n",
      " [27.49196964 37.65383533]\n",
      " [44.90882102 14.73518658]\n",
      " [15.36645123 12.54285363]\n",
      " [41.63046104 18.62457772]\n",
      " [36.05213095 38.8217706 ]\n",
      " [ 0.85219222 27.22190123]\n",
      " [19.97872788 12.14836792]\n",
      " [48.83986101 17.76409099]\n",
      " [52.97233893 11.47631846]\n",
      " [17.46830003 20.656496  ]\n",
      " [45.52581714 13.36554256]\n",
      " [27.65639033 47.14941755]\n",
      " [42.75312539 19.37832989]\n",
      " [44.49870547 19.25522105]\n",
      " [17.79951894 28.45571535]\n",
      " [ 2.94083873 30.488646  ]\n",
      " [43.21307182 16.96469176]\n",
      " [48.37588183 10.95393158]\n",
      " [35.0330945  39.50389997]\n",
      " [20.18659616 21.05096501]\n",
      " [45.98157432 17.86602315]\n",
      " [21.73498424 15.24527804]\n",
      " [20.85161455 23.02029469]\n",
      " [51.64591939 15.18873873]\n",
      " [46.31123218  8.38334668]\n",
      " [41.80897518 20.61786037]\n",
      " [46.86594544 10.53701837]\n",
      " [24.33335697 16.85377572]\n",
      " [26.22060704 11.17890643]\n",
      " [ 2.42181708 41.60985951]\n",
      " [40.06879371 12.09624594]\n",
      " [45.73560669 16.99363952]\n",
      " [37.29586391 45.16782437]\n",
      " [12.93830641 28.01536772]\n",
      " [33.4033128  46.88087069]\n",
      " [38.31886849 17.98774881]\n",
      " [20.2303988  22.41264684]\n",
      " [33.14415568 44.04670992]\n",
      " [22.48550431 18.20653748]\n",
      " [44.50106311 20.20260418]\n",
      " [50.92852634  4.66626357]\n",
      " [23.54499253 26.91989701]\n",
      " [37.4792409  -3.61953067]\n",
      " [25.6272157  42.32243291]\n",
      " [38.81325595 12.78543025]\n",
      " [22.24879131 21.92703094]\n",
      " [37.51286603 12.47000557]\n",
      " [51.14780396 18.63505294]\n",
      " [48.62593567 11.60276453]\n",
      " [50.87058176 17.03036299]\n",
      " [35.12087963 35.1159508 ]\n",
      " [51.61853264 20.43002412]\n",
      " [38.85008943 22.40051913]\n",
      " [27.96208852 36.42542967]\n",
      " [ 7.81282396 36.31393519]\n",
      " [32.60240161 43.74356762]\n",
      " [ 9.05045467 18.93936495]\n",
      " [35.6633267  33.76493449]\n",
      " [31.66345904 46.46911036]\n",
      " [39.39041856 19.97600514]\n",
      " [46.21521909 12.2187243 ]\n",
      " [31.02083995 38.97221492]\n",
      " [49.43838321 18.7974899 ]\n",
      " [26.31610451 41.51642869]\n",
      " [31.22499625 33.45636486]\n",
      " [23.1587153  12.32138786]\n",
      " [ 8.41036509 37.03707552]\n",
      " [43.92313626 22.78795159]\n",
      " [54.68803507 22.86500523]\n",
      " [-0.42768235 38.16823413]\n",
      " [40.88709203 15.54097659]\n",
      " [47.71037696 18.30262518]\n",
      " [46.27328318 16.10123533]\n",
      " [45.84349904  8.38798138]\n",
      " [48.59602654 10.19009774]\n",
      " [23.08913406 32.7446444 ]\n",
      " [31.69131396 11.75286217]\n",
      " [19.77345724 24.18068501]\n",
      " [20.21294552 23.76144   ]\n",
      " [43.98430703 13.21780412]\n",
      " [49.28277416  7.35027103]\n",
      " [38.74941573  8.68751586]\n",
      " [41.21114172 41.4187116 ]\n",
      " [30.45947233 40.26168941]\n",
      " [ 5.27480592 31.52593132]\n",
      " [11.12383689 23.21082845]\n",
      " [15.31002808 23.48522228]\n",
      " [19.31174301 12.06377332]\n",
      " [27.06829255 19.56522694]\n",
      " [55.73783166 22.42067354]\n",
      " [44.27889288 12.62513865]\n",
      " [46.97345266 19.00159824]\n",
      " [46.85623007 13.29325779]\n",
      " [53.94727683 10.77785835]\n",
      " [29.07263922 42.72236229]\n",
      " [37.54122814  5.62061282]\n",
      " [22.06236809 17.34849703]\n",
      " [ 9.9252419  22.72981498]\n",
      " [47.26721855  5.47999064]\n",
      " [22.51775728 20.23438409]\n",
      " [ 6.53830827 29.72425012]\n",
      " [30.7524091  47.88103645]\n",
      " [45.74899017  5.24366073]\n",
      " [44.9852174   7.93614261]\n",
      " [40.31942128 15.59224597]\n",
      " [30.46730585 42.2477504 ]\n",
      " [35.68137126  6.64906579]\n",
      " [34.34111648  6.53502565]\n",
      " [ 5.66805317 37.43694826]\n",
      " [44.04424599 12.28695626]\n",
      " [43.57917614 15.83566957]\n",
      " [12.37530315 22.39238715]\n",
      " [40.80498283 15.07270971]\n",
      " [50.64734694 20.05568548]\n",
      " [16.82570025 16.94750204]\n",
      " [33.72976339 12.57111739]\n",
      " [23.02755935 18.67825779]\n",
      " [38.30391476 16.16014688]\n",
      " [34.35166737 40.96275699]\n",
      " [29.98243574 43.92956005]\n",
      " [22.77455754 33.28203215]\n",
      " [24.58445747 25.74308856]\n",
      " [30.93877194 43.39342663]\n",
      " [ 4.57633254 27.00093975]\n",
      " [33.81655021 36.81770704]\n",
      " [24.42252681 17.7957153 ]\n",
      " [15.36674387 40.2651877 ]\n",
      " [46.07383423 20.18142021]\n",
      " [21.21547194  3.56808091]\n",
      " [56.41447749  7.38920864]\n",
      " [15.40328695  5.30955453]\n",
      " [35.50813819 30.95154917]\n",
      " [45.15761285  9.85646969]\n",
      " [43.47660799 10.50961992]\n",
      " [42.5824171  13.0353461 ]\n",
      " [40.49385384  8.17340793]\n",
      " [45.29994336 19.46119841]\n",
      " [17.02945854 17.67223195]\n",
      " [43.10845576 14.11397533]\n",
      " [34.37418262 26.08570212]\n",
      " [ 1.54836383 36.9765238 ]\n",
      " [-0.79834945 36.57511878]\n",
      " [52.55826748 19.78372241]\n",
      " [41.39425434 11.01442299]\n",
      " [55.06257446 11.38049018]\n",
      " [19.90764798 18.99493023]\n",
      " [26.97973003 48.74237911]\n",
      " [14.13839229 17.91087363]\n",
      " [15.03191434 10.02431859]\n",
      " [42.76655307 17.51259875]\n",
      " [34.51267055  9.53058987]\n",
      " [36.99962946 15.52451223]\n",
      " [43.11695306 16.26040704]\n",
      " [15.03506303 18.57337061]\n",
      " [42.87853867 14.31877688]\n",
      " [22.6628644  16.88488018]\n",
      " [-1.83349388 40.87535909]\n",
      " [34.60608755 34.7445121 ]\n",
      " [24.70241684 28.17866324]\n",
      " [44.90463895 10.54759544]\n",
      " [18.71774612 34.9159218 ]\n",
      " [23.78260869 19.02290786]\n",
      " [52.0780363  25.83795419]\n",
      " [36.03518243 35.43208465]\n",
      " [42.84316259  6.23175796]\n",
      " [43.37935126 19.35917965]\n",
      " [46.10093836  9.88961029]\n",
      " [ 9.61845129 35.06577   ]\n",
      " [46.9242286  18.48515145]\n",
      " [18.66376257 20.1152466 ]\n",
      " [43.58744778 18.8510222 ]\n",
      " [43.71743907  7.98208622]\n",
      " [16.91836261 31.53123399]\n",
      " [49.72925641 12.27310439]\n",
      " [19.28182814 48.59306077]\n",
      " [43.64853665 13.08196521]\n",
      " [40.00725943 21.64136839]\n",
      " [40.45152575 10.34476335]\n",
      " [52.0369917  15.04575629]\n",
      " [35.33322935 14.33133072]\n",
      " [18.27177086 10.54005447]\n",
      " [16.8821376  20.71245848]\n",
      " [17.92133024 18.66678649]\n",
      " [42.22795875  6.58138171]\n",
      " [43.33600538 12.87437941]\n",
      " [21.30184997 35.2965763 ]\n",
      " [49.36871032 14.80442115]\n",
      " [19.13645807 15.52112601]\n",
      " [40.36253813  5.87366868]\n",
      " [31.35734205 39.34149633]\n",
      " [23.25459316 45.34183792]\n",
      " [44.20424587  9.95323877]\n",
      " [48.77831674  6.77144794]\n",
      " [41.42534744 18.6122903 ]\n",
      " [48.43252012 11.2406488 ]\n",
      " [36.05037757 47.69036101]\n",
      " [47.99057536 16.28219515]\n",
      " [44.41697984 10.7259089 ]\n",
      " [25.0574604   9.1824057 ]\n",
      " [49.68518381  6.8487074 ]\n",
      " [23.82680385 28.92024038]\n",
      " [31.84945118 42.20065374]\n",
      " [29.92712336 34.45835568]\n",
      " [36.99597309 44.00943125]\n",
      " [37.37207499 16.9639059 ]\n",
      " [43.44824494 26.74288237]\n",
      " [45.46406667 14.35396634]\n",
      " [17.88939625 27.39713468]\n",
      " [44.25291282 22.14504224]\n",
      " [42.57963142  9.00015024]\n",
      " [47.57940956  9.21700022]\n",
      " [ 5.60626347 36.29804038]\n",
      " [ 7.15085636 32.44341538]\n",
      " [-1.15833808 36.29673242]\n",
      " [31.2984723  43.62889074]\n",
      " [44.28206723  8.20757003]\n",
      " [20.4190132  19.74903197]\n",
      " [16.62767682 26.38774859]\n",
      " [23.17248109 17.81401912]\n",
      " [44.39519134  8.31548948]\n",
      " [13.14849606 29.13746588]\n",
      " [36.22951497 13.66663121]\n",
      " [ 1.56296565 27.61207699]\n",
      " [49.28246885 21.72341306]\n",
      " [51.06326071 12.86537286]\n",
      " [25.66996226 34.65870067]\n",
      " [27.37320656 18.32371243]]\n",
      "初始化的mu1，mu2，mu3，mu4： [[0.02590709 0.44153678]\n",
      " [0.35933817 0.07102081]\n",
      " [0.32757448 0.50439989]\n",
      " [0.1791626  0.89968718]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "隐藏变量：\n",
      " [[0.12583345 0.09106637 0.25709003 0.52601015]\n",
      " [0.18666749 0.1583398  0.26912266 0.38587004]\n",
      " [0.15714761 0.15139742 0.2831503  0.40830468]\n",
      " ...\n",
      " [0.11570459 0.26253185 0.33912278 0.28264078]\n",
      " [0.13610925 0.10251515 0.26228304 0.49909256]\n",
      " [0.15876048 0.1859173  0.29564499 0.35967723]]\n",
      "迭代次数: 1\n",
      "估计的均值: [[30.14661682 22.16407259]\n",
      " [36.70677504 17.81370614]\n",
      " [34.23353566 20.97069867]\n",
      " [29.81381944 25.26147855]]\n",
      "估计的混合项系数: [0.14054908994497714, 0.18187661971728497, 0.29257799862569284, 0.3849962917120452]\n",
      "隐藏变量：\n",
      " [[1.70453020e-02 9.65653895e-06 2.87921370e-03 9.80065828e-01]\n",
      " [3.70162631e-01 1.25803527e-04 6.96865480e-03 6.22742911e-01]\n",
      " [2.30860814e-01 2.12843361e-03 5.12219769e-02 7.15788776e-01]\n",
      " ...\n",
      " [4.24271080e-05 9.69753044e-01 3.01967783e-02 7.75071399e-06]\n",
      " [3.33143056e-02 2.05847445e-05 4.14192506e-03 9.62523185e-01]\n",
      " [3.53648028e-01 5.25391877e-02 2.56499639e-01 3.37313146e-01]]\n",
      "迭代次数: 2\n",
      "估计的均值: [[20.10813865 18.46980378]\n",
      " [45.07766742 13.71926708]\n",
      " [39.56315002 18.5792746 ]\n",
      " [19.74270391 33.11016979]]\n",
      "估计的混合项系数: [0.08823179007256181, 0.4173921172617596, 0.1038226916597239, 0.3905534010059557]\n",
      "隐藏变量：\n",
      " [[1.46910598e-06 2.28646367e-12 1.49183368e-07 9.99998382e-01]\n",
      " [9.87206668e-01 1.00033266e-13 2.26527316e-09 1.27933299e-02]\n",
      " [6.52310929e-01 4.61785835e-07 9.03512225e-04 3.46785097e-01]\n",
      " ...\n",
      " [3.36080920e-15 9.96562269e-01 3.43773053e-03 2.31492245e-20]\n",
      " [4.52289265e-05 5.91693382e-12 2.69633743e-07 9.99954501e-01]\n",
      " [9.51372668e-01 3.74219809e-04 4.58499345e-02 2.40317789e-03]]\n",
      "迭代次数: 3\n",
      "估计的均值: [[19.38655976 17.83099339]\n",
      " [45.66206882 13.62092279]\n",
      " [39.76068466 20.22426526]\n",
      " [19.23032053 36.61034903]]\n",
      "估计的混合项系数: [0.16612881475244926, 0.41092824841373515, 0.11151992166751651, 0.311423015166299]\n",
      "隐藏变量：\n",
      " [[6.06622789e-07 7.36369012e-13 8.22059114e-07 9.99998571e-01]\n",
      " [9.99844354e-01 1.25358624e-14 7.96091673e-10 1.55645606e-04]\n",
      " [9.76314520e-01 2.41432738e-07 1.63080811e-03 2.20544302e-02]\n",
      " ...\n",
      " [1.42454577e-15 9.98301215e-01 1.69878490e-03 3.00401274e-23]\n",
      " [5.14822058e-05 3.82573451e-12 2.08421550e-06 9.99946434e-01]\n",
      " [9.70518411e-01 1.39620355e-04 2.93176467e-02 2.43215292e-05]]\n",
      "迭代次数: 4\n",
      "估计的均值: [[19.09767172 18.62730559]\n",
      " [45.63152603 13.4102973 ]\n",
      " [39.80751009 21.64437486]\n",
      " [18.99328489 37.5880376 ]]\n",
      "估计的混合项系数: [0.18849120745598896, 0.4076001578394058, 0.12187149826304906, 0.28203713644155604]\n",
      "隐藏变量：\n",
      " [[1.98587738e-06 6.30533624e-13 5.23863837e-06 9.99992775e-01]\n",
      " [9.99963732e-01 9.03294206e-15 5.30091988e-10 3.62676355e-05]\n",
      " [9.92670983e-01 1.60632935e-07 1.64372515e-03 5.68513157e-03]\n",
      " ...\n",
      " [6.65888067e-16 9.99094445e-01 9.05554960e-04 3.42052635e-24]\n",
      " [1.77259466e-04 4.21686667e-12 1.14744207e-05 9.99811266e-01]\n",
      " [9.75051714e-01 1.38440416e-04 2.48039872e-02 5.85840658e-06]]\n",
      "迭代次数: 5\n",
      "估计的均值: [[18.87349598 19.05905892]\n",
      " [45.49117023 13.32031841]\n",
      " [39.85482458 22.9200505 ]\n",
      " [18.78303055 38.01058141]]\n",
      "估计的混合项系数: [0.19831824645657986, 0.41194621110574986, 0.12358153038180082, 0.2661540120558687]\n",
      "隐藏变量：\n",
      " [[3.79567875e-06 7.69169107e-13 2.35548734e-05 9.99972649e-01]\n",
      " [9.99981046e-01 1.01936783e-14 3.44161748e-10 1.89538746e-05]\n",
      " [9.95381126e-01 1.67116555e-07 1.61857462e-03 3.00013183e-03]\n",
      " ...\n",
      " [3.79485706e-16 9.99556191e-01 4.43809078e-04 1.05843761e-24]\n",
      " [3.38684052e-04 5.75112659e-12 4.05788661e-05 9.99620737e-01]\n",
      " [9.80820973e-01 1.77128148e-04 1.89987796e-02 3.11972023e-06]]\n",
      "迭代次数: 6\n",
      "估计的均值: [[18.76825601 19.27171007]\n",
      " [45.33290539 13.35865502]\n",
      " [39.75717102 24.14233074]\n",
      " [18.57770093 38.21640352]]\n",
      "估计的混合项系数: [0.20316859104477053, 0.42197582593257105, 0.11808169636320132, 0.2567738866594569]\n",
      "隐藏变量：\n",
      " [[5.61018007e-06 1.18562692e-12 9.45006857e-05 9.99899889e-01]\n",
      " [9.99985852e-01 1.38315128e-14 2.66215040e-10 1.41475832e-05]\n",
      " [9.96254825e-01 2.08795586e-07 1.65554163e-03 2.08942444e-03]\n",
      " ...\n",
      " [2.94273249e-16 9.99828509e-01 1.71490802e-04 4.81066501e-25]\n",
      " [4.92525028e-04 9.16103769e-12 1.25864248e-04 9.99381611e-01]\n",
      " [9.86182281e-01 2.34444351e-04 1.35811266e-02 2.14802627e-06]]\n",
      "迭代次数: 7\n",
      "估计的均值: [[18.74031121 19.36419911]\n",
      " [45.1941791  13.49386578]\n",
      " [39.36956954 25.57937772]\n",
      " [18.29369121 38.30561058]]\n",
      "估计的混合项系数: [0.20533275825395486, 0.435698906741599, 0.10919890822447552, 0.24976942677997035]\n",
      "隐藏变量：\n",
      " [[7.69263256e-06 2.20871690e-12 5.22753109e-04 9.99469554e-01]\n",
      " [9.99986967e-01 1.96889517e-14 2.81278825e-10 1.30323195e-05]\n",
      " [9.96415283e-01 2.78819406e-07 1.96958431e-03 1.61485389e-03]\n",
      " ...\n",
      " [2.76428158e-16 9.99961649e-01 3.83511622e-05 2.22821424e-25]\n",
      " [6.50249752e-04 1.65282667e-11 5.12950288e-04 9.98836800e-01]\n",
      " [9.90348660e-01 3.03272682e-04 9.34647513e-03 1.59186570e-06]]\n",
      "迭代次数: 8\n",
      "估计的均值: [[18.76949962 19.385985  ]\n",
      " [45.06649096 13.70336083]\n",
      " [38.54119405 27.8070144 ]\n",
      " [17.72370679 38.29335346]]\n",
      "估计的混合项系数: [0.2059613187002078, 0.45253463120838844, 0.1006229311124549, 0.2408811189789486]\n",
      "隐藏变量：\n",
      " [[1.20742163e-05 5.60790990e-12 7.38008868e-03 9.92607837e-01]\n",
      " [9.99985270e-01 2.90713779e-14 3.57158124e-10 1.47298040e-05]\n",
      " [9.96268287e-01 3.89101803e-07 2.53015610e-03 1.20116756e-03]\n",
      " ...\n",
      " [2.99253036e-16 9.99997623e-01 2.37654146e-06 6.39274157e-26]\n",
      " [9.30981319e-04 3.74477588e-11 4.42576404e-03 9.94643255e-01]\n",
      " [9.94965750e-01 3.86081049e-04 4.64709032e-03 1.07842861e-06]]\n",
      "迭代次数: 9\n",
      "估计的均值: [[18.88692163 19.34437711]\n",
      " [44.9143991  13.99257451]\n",
      " [36.97515853 31.49307494]\n",
      " [16.55515619 38.13680956]]\n",
      "估计的混合项系数: [0.2053268430906146, 0.47328688559376486, 0.09556242852776323, 0.22582384278785728]\n",
      "隐藏变量：\n",
      " [[2.16264230e-05 1.97011095e-11 2.86513543e-01 7.13464831e-01]\n",
      " [9.99978377e-01 4.82457346e-14 3.49400730e-10 2.16222858e-05]\n",
      " [9.97134663e-01 5.88071825e-07 2.15366929e-03 7.11079350e-04]\n",
      " ...\n",
      " [4.05150798e-16 9.99999991e-01 9.42893353e-09 5.75840633e-27]\n",
      " [1.70738010e-03 1.29458298e-10 9.24918011e-02 9.05800819e-01]\n",
      " [9.98747498e-01 4.94184612e-04 7.57790670e-04 5.26450132e-07]]\n",
      "迭代次数: 10\n",
      "估计的均值: [[19.14778629 19.23278211]\n",
      " [44.74399901 14.32913915]\n",
      " [34.17138437 36.80375221]\n",
      " [13.39045192 37.31828154]]\n",
      "估计的混合项系数: [0.2032198519959385, 0.4945287620525424, 0.11268433752065023, 0.18956704843086897]\n",
      "隐藏变量：\n",
      " [[3.54547885e-06 7.00853147e-12 9.93676213e-01 6.32024122e-03]\n",
      " [9.99939306e-01 8.85895010e-14 1.41822705e-10 6.06942440e-05]\n",
      " [9.99370308e-01 9.10286381e-07 4.69398228e-04 1.59383722e-04]\n",
      " ...\n",
      " [7.76676434e-16 1.00000000e+00 3.47084333e-13 9.79335502e-30]\n",
      " [1.82403661e-03 2.87118345e-10 8.96738791e-01 1.01437172e-01]\n",
      " [9.99386644e-01 6.00028560e-04 1.32468641e-05 8.08181933e-08]]\n",
      "迭代次数: 11\n",
      "估计的均值: [[19.4785959  19.13594143]\n",
      " [44.65098158 14.55726618]\n",
      " [31.59052605 39.37169723]\n",
      " [ 8.53430927 35.73074171]]\n",
      "估计的混合项系数: [0.20077642817736796, 0.5062218625605244, 0.14714729011396438, 0.1458544191481434]\n",
      "隐藏变量：\n",
      " [[1.18142402e-06 3.61407679e-12 9.99993372e-01 5.44632467e-06]\n",
      " [9.99895478e-01 1.34201940e-13 2.16810982e-10 1.04522263e-04]\n",
      " [9.99791177e-01 1.12800688e-06 2.02378435e-04 5.31668681e-06]\n",
      " ...\n",
      " [1.67579027e-15 1.00000000e+00 2.75309895e-16 2.20698626e-34]\n",
      " [8.32377396e-04 1.99535035e-10 9.98806816e-01 3.60806868e-04]\n",
      " [9.99387427e-01 6.11295035e-04 1.27572083e-06 1.77227828e-09]]\n",
      "迭代次数: 12\n",
      "估计的均值: [[19.61665311 19.24500078]\n",
      " [44.62203559 14.59561056]\n",
      " [30.62796198 39.50238004]\n",
      " [ 6.54026083 35.37178394]]\n",
      "估计的混合项系数: [0.2025969621372135, 0.5079306668076163, 0.16086938065460454, 0.12860299040056536]\n",
      "隐藏变量：\n",
      " [[1.10713423e-06 2.98152704e-12 9.99998666e-01 2.27190910e-07]\n",
      " [9.99943692e-01 1.49807713e-13 6.05232825e-10 5.63078824e-05]\n",
      " [9.99730272e-01 1.10208031e-06 2.68037399e-04 5.88590317e-07]\n",
      " ...\n",
      " [2.18757602e-15 1.00000000e+00 6.72553605e-17 1.04186945e-36]\n",
      " [6.69182759e-04 1.47117028e-10 9.99311183e-01 1.96339404e-05]\n",
      " [9.99404495e-01 5.94136723e-04 1.36853247e-06 1.56766017e-10]]\n",
      "迭代次数: 13\n",
      "估计的均值: [[19.59059389 19.33559656]\n",
      " [44.62068879 14.59908462]\n",
      " [30.30532933 39.51726388]\n",
      " [ 5.93073742 35.3091522 ]]\n",
      "估计的混合项系数: [0.20452134148870793, 0.5080294012741066, 0.16522412324106667, 0.12222513399611841]\n",
      "隐藏变量：\n",
      " [[1.13798854e-06 2.76651880e-12 9.99998780e-01 8.17067930e-08]\n",
      " [9.99959788e-01 1.47316477e-13 8.43943940e-10 4.02107123e-05]\n",
      " [9.99699607e-01 1.07644820e-06 2.99046036e-04 2.70749488e-07]\n",
      " ...\n",
      " [2.01289501e-15 1.00000000e+00 4.32304251e-17 1.76084503e-37]\n",
      " [6.47267207e-04 1.30678093e-10 9.99345236e-01 7.49690355e-06]\n",
      " [9.99396916e-01 6.01598379e-04 1.48562850e-06 6.84352228e-11]]\n",
      "迭代次数: 14\n",
      "估计的均值: [[19.56516591 19.37992018]\n",
      " [44.61917735 14.6002512 ]\n",
      " [30.20565972 39.51402559]\n",
      " [ 5.73168343 35.31980372]]\n",
      "估计的混合项系数: [0.2054174593975447, 0.5081035552876539, 0.16651061404429524, 0.11996837127050615]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "隐藏变量：\n",
      " [[1.16331181e-06 2.71055026e-12 9.99998778e-01 5.85831137e-08]\n",
      " [9.99965636e-01 1.45770841e-13 9.35955133e-10 3.43633990e-05]\n",
      " [9.99687505e-01 1.07009176e-06 3.11220943e-04 2.03647748e-07]\n",
      " ...\n",
      " [1.88141803e-15 1.00000000e+00 3.81845957e-17 9.33898936e-38]\n",
      " [6.47257940e-04 1.26100671e-10 9.99347307e-01 5.43495071e-06]\n",
      " [9.99388330e-01 6.10118412e-04 1.55148181e-06 5.04613386e-11]]\n",
      "迭代次数: 15\n",
      "估计的均值: [[19.55237788 19.39960906]\n",
      " [44.61833299 14.60066759]\n",
      " [30.17443186 39.51192323]\n",
      " [ 5.66133125 35.33116983]]\n",
      "估计的混合项系数: [0.20580431411168812, 0.5081414474635888, 0.16691226610537085, 0.11914197231935207]\n",
      "隐藏变量：\n",
      " [[1.17742128e-06 2.69483166e-12 9.99998770e-01 5.21645306e-08]\n",
      " [9.99967834e-01 1.45078433e-13 9.65861581e-10 3.21649757e-05]\n",
      " [9.99683469e-01 1.06813336e-06 3.15280343e-04 1.82860206e-07]\n",
      " ...\n",
      " [1.81997188e-15 1.00000000e+00 3.67985193e-17 7.37301328e-38]\n",
      " [6.49989563e-04 1.24745621e-10 9.99345159e-01 4.85132827e-06]\n",
      " [9.99383818e-01 6.14603367e-04 1.57874339e-06 4.49666329e-11]]\n",
      "迭代次数: 16\n",
      "估计的均值: [[19.54692984 19.40791682]\n",
      " [44.61795719 14.60081775]\n",
      " [30.16422264 39.51125273]\n",
      " [ 5.63564528 35.33690981]]\n",
      "估计的混合项系数: [0.2059652655141371, 0.508157790338252, 0.16704274170356428, 0.11883420244404583]\n",
      "隐藏变量：\n",
      " [[1.18423308e-06 2.69015195e-12 9.99998766e-01 5.00239239e-08]\n",
      " [9.99968670e-01 1.44792537e-13 9.75188736e-10 3.13294162e-05]\n",
      " [9.99682236e-01 1.06736398e-06 3.16521036e-04 1.75543783e-07]\n",
      " ...\n",
      " [1.79442096e-15 1.00000000e+00 3.63558833e-17 6.74647708e-38]\n",
      " [6.51906880e-04 1.24327836e-10 9.99343437e-01 4.65555178e-06]\n",
      " [9.99381864e-01 6.16547026e-04 1.58870113e-06 4.30464406e-11]]\n",
      "迭代次数: 17\n",
      "估计的均值: [[19.54475755 19.41130051]\n",
      " [44.61780271 14.60087408]\n",
      " [30.16075642 39.51108912]\n",
      " [ 5.62612861 35.33939031]]\n",
      "估计的混合项系数: [0.2060303869932675, 0.5081644320724562, 0.16708637227082568, 0.1187188086634507]\n",
      "隐藏变量：\n",
      " [[1.18725303e-06 2.68870445e-12 9.99998763e-01 4.92597054e-08]\n",
      " [9.99968988e-01 1.44679360e-13 9.78120744e-10 3.10105643e-05]\n",
      " [9.99681872e-01 1.06703743e-06 3.16887699e-04 1.72845834e-07]\n",
      " ...\n",
      " [1.78428494e-15 1.00000000e+00 3.62028519e-17 6.52452518e-38]\n",
      " [6.52895905e-04 1.24194925e-10 9.99342518e-01 4.58554850e-06]\n",
      " [9.99381080e-01 6.17328228e-04 1.59219760e-06 4.23412465e-11]]\n",
      "迭代次数: 18\n",
      "估计的均值: [[19.54391538 19.41264595]\n",
      " [44.61774163 14.60089549]\n",
      " [30.15954463 39.51105995]\n",
      " [ 5.62257555 35.34040001]]\n",
      "估计的混合项系数: [0.2060562090864686, 0.5081670471897993, 0.16710134482356645, 0.11867539890016598]\n",
      "隐藏变量：\n",
      " [[1.18852085e-06 2.68824022e-12 9.99998762e-01 4.89791015e-08]\n",
      " [9.99969110e-01 1.44635400e-13 9.79064740e-10 3.08889204e-05]\n",
      " [9.99681766e-01 1.06689941e-06 3.16995441e-04 1.71833840e-07]\n",
      " ...\n",
      " [1.78034976e-15 1.00000000e+00 3.61477323e-17 6.44271519e-38]\n",
      " [6.53344441e-04 1.24151306e-10 9.99342096e-01 4.55984090e-06]\n",
      " [9.99380774e-01 6.17632155e-04 1.59342004e-06 4.20773855e-11]]\n",
      "迭代次数: 19\n",
      "估计的均值: [[19.54359295 19.41317235]\n",
      " [44.61771796 14.60090365]\n",
      " [30.1591116  39.51105893]\n",
      " [ 5.62124309 35.34079922]]\n",
      "估计的混合项系数: [0.20606630141048066, 0.5081680593841258, 0.1671066011293074, 0.11865903807608584]\n",
      "隐藏变量：\n",
      " [[1.18903472e-06 2.68808619e-12 9.99998762e-01 4.88747389e-08]\n",
      " [9.99969156e-01 1.44618492e-13 9.79377444e-10 3.08425914e-05]\n",
      " [9.99681734e-01 1.06684250e-06 3.17027449e-04 1.71451935e-07]\n",
      " ...\n",
      " [1.77883852e-15 1.00000000e+00 3.61274566e-17 6.41209832e-38]\n",
      " [6.53534431e-04 1.24136562e-10 9.99341915e-01 4.55028223e-06]\n",
      " [9.99380657e-01 6.17748687e-04 1.59385154e-06 4.19779683e-11]]\n",
      "迭代次数: 20\n",
      "估计的均值: [[19.54347022 19.41337609]\n",
      " [44.61770888 14.60090675]\n",
      " [30.15895433 39.51106139]\n",
      " [ 5.620742   35.34095453]]\n",
      "估计的混合项系数: [0.20607020603560638, 0.5081684475104998, 0.16710848166181652, 0.11865286479207769]\n",
      "隐藏变量：\n",
      " [[1.18923830e-06 2.68803353e-12 9.99998762e-01 4.88356705e-08]\n",
      " [9.99969174e-01 1.44612026e-13 9.79483952e-10 3.08249787e-05]\n",
      " [9.99681725e-01 1.06681959e-06 3.17037165e-04 1.71307524e-07]\n",
      " ...\n",
      " [1.77826159e-15 1.00000000e+00 3.61199145e-17 6.40057162e-38]\n",
      " [6.53611745e-04 1.24131446e-10 9.99341841e-01 4.54670513e-06]\n",
      " [9.99380613e-01 6.17793069e-04 1.59400593e-06 4.19404149e-11]]\n",
      "迭代次数: 21\n",
      "估计的均值: [[19.54342364 19.41345436]\n",
      " [44.61770541 14.60090793]\n",
      " [30.15889654 39.51106309]\n",
      " [ 5.62055323 35.34101436]]\n",
      "估计的混合项系数: [0.20607170605691397, 0.5081685955554439, 0.167109164598591, 0.11865053378905156]\n",
      "0.00043761450425083837 4.662006052258416e-06\n"
     ]
    }
   ],
   "source": [
    "if __name__ == '__main__':\n",
    "    iter_num=1000  #迭代次数\n",
    "    N=500         #样本数目\n",
    "    k=4            #高斯模型数\n",
    "    probility = np.zeros(N)    #混合高斯分布\n",
    "    u1=[5,35]\n",
    "    u2=[30,40]\n",
    "    u3=[20,20]\n",
    "    u4=[45,15]\n",
    "    sigma=np.matrix([[30, 0], [0, 30]])               #协方差矩阵\n",
    "    alpha=[0.1,0.2,0.3,0.4]         #混合项系数\n",
    "    generate_data(sigma,N,u1,u2,u3,u4,alpha)     #生成数据\n",
    "    #迭代计算\n",
    "    for i in range(iter_num):\n",
    "        err=0     #均值误差\n",
    "        err_alpha=0    #混合项系数误差\n",
    "        Old_mu = copy.deepcopy(mu)\n",
    "        Old_alpha = copy.deepcopy(alpha_)\n",
    "        \n",
    "        e_step(sigma,k,N)     # E步\n",
    "        m_step(k,N)           # M步\n",
    "        \n",
    "        print(\"迭代次数:\",i+1)\n",
    "        print(\"估计的均值:\",mu)\n",
    "        print(\"估计的混合项系数:\",alpha_)\n",
    "        for z in range(k):\n",
    "            err += (abs(Old_mu[z,0]-mu[z,0])+abs(Old_mu[z,1]-mu[z,1]))      #计算误差\n",
    "            err_alpha += abs(Old_alpha[z]-alpha_[z])\n",
    "        if (err<=0.001) and (err_alpha<0.001):     #达到精度退出迭代\n",
    "            print(err,err_alpha)\n",
    "            break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD3CAYAAADi8sSvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eXhcV5kn/HtrL1WVVCWVLGuX5TWJl9iWt5AMaUw6JEBiNxB6oTuBDoHmm2+mv6ZZpme+7maaXufrNHT3PB0CfAQGpoFAnEAgIYkhiU2cSPImyfGqfZdKKtW+3zN/vPfUvSqXpJItR5Zd7/PUU3Xrbufe+97fec/vXQ4JIVCUohSlKEW5McSw3A0oSlGKUpSiLJ0UQb0oRSlKUW4gKYJ6UYpSlKLcQFIE9aIUpShFuYGkCOpFKUpRinIDSRHUi1KUohTlBpIbHtSJ6C+J6LvL3Y4bSYjoESI6uojt+4jovdeyTdeLLPbeXMHxXyCih3XLXyYiHxGNEVEDEYWJyHgFx20iIkFEpgK3f4qIvrzY8yyVXE86RUSvEtGjy90OKTc8qBflcrnelFAvKrCsW+52XK8ihLhPCPFtACCiegCfBXCrEGK1EGJACOEUQmSWt5Wz5Wr1bbk7kKUUIrqbiBS189V/9qnrX1XfgW05+z2r/n/3Que4rkC9UCuhKHNL8R7eVNIIYEoIMbHcDbme5Tp8J0bUzlf/OaZbfwHAH8gFIqoAsBfAZCEHX3ZQV4dRXyCiDgARIjIR0ReJqJuIQkT0NhEd1G3/CBEdJaL/j4j8RNRLRPfp1q8hotfUfV8G4M053wNEdIaIZtRe8ZactnyOiDqIKEJE3ySiKnXIGyKiV4jIM8+1fJ6IRolohIge1VudRGRV2zxARONE9AQR2dV1dxPREBF9logm1GN8XHfcQvb9AhGNAfgWEXmI6HkimlTv0fNEVKdu/9cA7gLwr6qF8K/q/5uI6GUimiai80T0kO78FUT0EyIKElErgLULPNPfJ6J+Ipoiov+as243ER1T7/8oEf0rEVnUda+rm51W2/bR+a5lOYWI6onoGbVdU/I+5tnuq0Q0qN6740R0l27dbiJqV9eNE9Hj6v82IvquetwZImojoip13auqbr0XwMsAatR79RTlUChEVKbq8CgRDRNTNUZ1nVHVKR8R9QB4/wLXu52ITqjvwQ8A2HTrrkTf5rwvOed9DMDvAfi8uv9PdatvJ35XA0T0AyKyqftc9k6o/3+SiC6pOv4TIqpR/7+MeiLd6EK9V/+o3qteIvqPudsDaCSiX6v35yUimoU7i5TvAfgoaTTa7wA4BCBZ0N5CiGX9AOgDcApAPQC7+t9HANSAO52PAogAqFbXPQIgBeCTAIwA/gjACABS1x8D8DgAK4D/ACAE4Lvqug3qse4BYAbweQCXAFh0bXkTQBWAWgATAE4A2K4e75cA/mKO63gfgDEAtwEoAfC/AAgA69T1XwHwEwDlAFwAfgrgb9V1dwNIA/jvarvuBxAF4FnEvn+vttEOoALAh9R2uAA8DeBZXVtfBfCobtkBYBDAxwGYAOwA4ANwm7r++wB+qG63GcAwgKNz3IdbAYTVe29Vn0UawHvV9TvBVocJQBOAswD+WLd/9p6py/NeyzLprBHAaQD/pN4TG4A7dfp5VLftx9RrMIGpkjEANp2u/r762wlgr/r7U+ozLlHPtRNAae6zU5/9kO5cTer9M6nLzwL4mtrGVQBaAXxKXfdpAOfA7105gF/p9825XguAfgD/D1g/Pwx+B79cyDPK1beF7kue8z8lz5WDG61gnChX9ejT87wT7wHr9A71v38B8Hq++5bnPn8awNsA6gB4ALySc59fBdANxhe7uvx3c1zLrGeWZ/2rAB4F8BKA+9T/WgHsAzAE4O4F9XM5Xw7dw/nEAtucAvCg7qW5pFtXot7g1QAa1Ifp0K3/39BA/f8F8EPdOgMYoO7WteX3dOt/DODfdMv/N+YAFAD/P1SgVZfXqe1aB4DAncla3fp9AHp1DzqWo1QTYPArZN8k5ngh1G1uB+Cf6yUDd5xHcvb5GoC/AINKCsAm3bq/wdyg/ucAvq9bdqjte+8c2/8xgEO65VmgvtC1LJPO7gMPhfMB4CNz3Rt1vR/ANvX36wC+BMCbs80nALwBYGue/bPPDvOAOtgwSUA1lNT1vwPgV+rvX0IFQXX5NzE3qP8H6Awn9b83kAO0herbQvclz7qncs8Fflc/plv+BwBPzPVOAPgmgH/QLTtVvW7CwqD+S6idobr8XlwO6v9Nt/4zAF6c41ruBqAAmMn5OPTnBXd6/w5gI4AL6rqCQP164ZoG9QtE9AcA/gR8swF+APrhzJj8IYSIEpF+G78QIqLbth9sjQDcq/fr9lWIaBBslUsZ1/2O5Vl2znENNQDa57imSnDnc1xtK8BgrY9SmBJCpHXLUfVchew7KYSIZ1cSlYCtyPeBLQsAcBGRUeR3ojUC2ENEM7r/TODRRqX6W389/ZhbavTbCiEiRDSla9sGsPXeol6XCcDxuQ52BdfyTkg9gP6c55VXiOiz4Je0BgwEpdB0+Q/Bo7NzRNQL4EtCiOfB970ewPeJyA3guwD+qxAitYg2NoKt6lGd3higPZtZzwkLP9NhoSJL7vZX8owWuC+Fypjud1Q9lpRZ74S67oRcEEKEVb2sBRt280nuvRrMs01uW+bCCYA59YUoxGcA/COAKbA+FCzLzqmrklUWImoE8HUA/xFAhRDCDaALDGQLySgADxE5dP816H6PgJVdnovAL89CD7UQGQUPz6TU6377wB3CbUIIt/opE0LM9+AXs6/I2eez4B5+jxCiFGxpAdo9zN1+EMBruuO7BTtv/ghskaZzrqcBc8uoflv1ha/Qrf838LB/vdq2P8P8z3aha1kOGQTQQAs44FSe+AsAHgJTaW4AAahtF0JcFEL8Dpga+XsAPyIihxAiJYT4khDiVgB3APgAdI6zRbQxAR4FyGdaKoS4TV0/6zlh4WdaS7reIWf7RenbQvclj+TqayGSu0/uu+8A6+UweCQMsJEhZbXu93zv9jURIUQUwAtgenlFgrpeHOAHMgkAxA7DzYXsKIToB1vLXyIiCxHdCeCDuk1+COD9RLSfiMxgZUyAh5JXKz8E8HEiukUFsj/XtUsBd1T/RESr1OuqJaJ7C7imK9nXBe4IZoioHEyj6GUcQLNu+XkAG4gdnGb1s4uIblEtrWcA/CURlRDRrQAenufcPwLwASK6k9gB+t8xW89cAIIAwkS0Cay087VtoWtZDmkFv+h/R0QOYsfmu/Js5wJ3iJMATET052CLFABARB8jokr1GctRUoaIfoOItqiOsiCYJljUqEQIMQrmZf+RiEqJyEBEa4no3eomPwTwn4iojtj5/8V5DndMvY7/RBzI8FsAdudc52L0bd77kkdy978S+d/g9/N2IrKCKcS3hBB9QohJMLh/THWKfgKzgwF+COA/q++dG9whvRPyZwDeLYToW8xO1x2oCyHeBg87joEf5hYAv17EIX4XwB4A02Dl+o7u2OfBXNW/gC3gDwL4oBCiMK/y/O1+AcA/gx1Ol9T2A9xpAKwIlwC8SURBsLNlY4GHX+y+XwE7bHxgx++LOeu/CuDDxJEK/yyECIE51d8GWzRj0JxMAI+anOr/T0GNJsgnQogzAP4v8Es0CuZKh3Sb/Cn4GYXAndUPcg7xlwC+TRz18VAB1/KOi9rRfRDsLxkAX99H82z6C7C1dQFMV8Qxe+j+PgBniCgMfia/rVIGq8GdYxDsAHwNTMEsVv4A7OR8G/wcfgSgWl33dbV9p8G0xDNzHUR9P34L7C/wg69Vv/2i9A0L35dc+SaAW1WdeHbeK577Gg6DfWo/BuvlWrC+S/kkgM+B6Y7bMNvQ+zq4g+wAcBLAz8Gd0pXSfzJiSf/5UJ42jwghFp3IJiNGirLEQhwq2QXAWgj3WpSiFGVlCHEI9RNCiMYFN14Gue4s9ZUsRHRQpX08YEv3p0VAL0pRVrYQkZ2I7lepp1owA3Bouds1lxRBfWnlU2CesBs8NMvli4tSlKKsPCFw6KkfTL+chc5ndr1JkX4pSlGKUpQbSIqWelGKUpSi3EDyjiYfeb1e0dTU9E6esig3kRw/ftwnhKh8p89b1OuiXGtZjG6/o6De1NSE9vb2hTcsSlGuQIhovqzIayZFvS7KtZbF6Pb1UiZgWUQIYHAQ8PkArxeorwdoOfMUi1KUohTlKuWmBXUhgEOHgLY2wGAAFAXYtQs4eLAI7EUpSlFWrty0oD44yIDe2KiBelsb0NICNMxXBaMoRVkhIoRAIjGIVMoHs9kLq7UeVLRYbni5aUHd52MwN6jxP/K3z7d0oL4U9E6RIirKlYgQAj7fIQSDbeAgNwWlpbvg9R4sAvsNLjctqHu9bJ0rimapKwr/vxSyFPROIccogn5R8kkiMYhgsA02WyOIDBBCQTDYBperBTZbcSh6I8tNC+r19QyQuYBZv0RFNZeC3lnoGFfScRQ7gZtDUikfAAOIeCjK3wakUr4lAXUhBAYHB+Hz+eD1elFfv3hqZymOUZTL5aYFdSLgwAGguhro6QGam4E9e5YO4JaC3tEfQwhgZgYYHwc6O3m5qws4fBjYvBkwGmeDfn395eANaJ1AKAREo8DevcCjj2rtLMqNIWazF4ACIZSspQ4o6v9XJ0IIHDp0CG1tbTAYDFAUBbt27cLBg4VTO4Ucowj6VyY3LagLATz7rGbldnUBo6NXF/2it4ITiaundyRFlMkAp08DfX1AMAj8y78AJSWA2w2cO8frb799Nui3tQHt7bMt+J07+f/paWBggM/xrW9xux97rGix30hitdajtHTXZZy61Xr1Q9HBwUG0tbWhsbExC8htbW1oaWlBQ4EWizxGQ0MDAoEAQqEQDh8+jJ07d6KxsXHRHUexA9DkpgX1pY5+yaVCMhk+Zn//ldM7kiI6fBjo6ABcLqCqChgbA2IxYP16YGQE6O0FwmG+Jr8fCAT4fC0tvJ3Lxe2yWtlCHxgAyssZxIUA3nwTuO++YtTPjSREhIqKA7BYqhGP98Bma4bLtWdJgM7n88FgMMCgDu/kb5/PVzCo+3w+EBFOnz6N/v5+EBFmZmbw9NNP4yMf+Qi6urpw+PBhbN68GQaDAVNTU3j++edRXV2NPXv2YGhoKAvgdXV1ePbZZ9HW1oZQKIRoNIq9e/fi0UcfzbbxZpKbFtSXOvolXyfR1wc88ACD6Vz89UIc986d2nqTiamiiQner6ICaGri8547xxSM2w0MDzNw+/1s0QvB33feyZQLoAF6PM7/dXYW+fUbSYQQmJp6Nmuph8NdSCZHryr6RVrDo6OjmJmZQV1dHYxGIxRFgaIo8C5iGOr1ehEIBNDX14eKCp7tMJPJ4JVXXkFPTw8SiQTOnj2LTIbnoTh37hz8fj9CoRCampoQiUQQj8dht9uxdetW9Pb2wu/3Y0Adgn7rW9+CEAKPPfbYTWex37SgvtTRL/k6CaORAX3Hjvz7zOfoBGbz35cusfXf1MSWOBHTKBL0o1Fg40bA4eC2BIN8/PJyPm5vL+DxMIcuKZdxdUrteBz4wQ+As2eBj3yEO7Wb7D244WSpo1/0dAgRYXR0FKOjo9i0aROEENi1axfqc4ah81EidXV1cLvdGBkZQSKRgNVqRUlJCbq7u7F27VqsW7cOw8PDOHv2LHw+H4QQiMfjGBgYwK9//WusX78eDocDiqLg7NmzqKurw+TkJMrLy0FEEELgzTffxH333Vfw6OFGkZsW1Jc6+uVKOon5KCBAW0fElMuxY2yJW9VJ5oJB4Fe/YovdZmOglvSL3CeVAsxmoKaG93v0Ud7u5ZfZonc4AIuFAf6ZZ3gk8N73FjNrV7osdfRLLo9eX1+Pzs5OvOtd78KWLVsu47Dn48QB4Nlnn8X09DTMZjMSiQQymQxmZmYwPT2NM2fOIJ1Oo6mpCa+99hqmpqZQWVmZ7TSCwSCICOXl5RBC4Ny5czAYDDCbzVlAj8fjiEQiOHz4MLZu3Zrd/2aw2m9aUCdi4GppWZrwvivpJOajgPTLAFvkvb0cpbNuHdMlb7zBoO3z8bq+Pt4+meRPIMD7VlVpHYzBAHzyk0zPXLjA64eGgLo6oKyMO41iZu3Kl6WOfsnl0Y1GIzweD6qrq/NawvM5UwGgra0NW7duhRACZ8+exaVLl1BXV4fKykqsXr0a58+fx2233Ybq6moIIbKW+fDwMIQQCIVCiEQisNvtsFgs2LRpE44fPw5FUTA0NIRUKoX+/n5cunQJzc3NWLduHXbv3r2oCJ2VKtc9qF9JXLXcZ3KSwc1iASorL9+XiIHrasBL376dO/kzNVVYWxey7vXr3G622m02trA7O/l3czM7S1Mp7ZpKSvjb5eL/iWa3Y2iII2VqavjYoRC3mQgoLeXlpcysLco7L0sd/eL1erPcuQTp+Xj0+ZypctloNGL79u0wm82YmZnBrl27kEwm0dbWhtHRUQwNDaGiogKKomBgYABOpxMjIyPIZDIYGRnB9PQ03G43mpqa8NnPfha/+MUv8J3vfAdTU1OwWCxIJBKoq6tDKBRCaWnpoiN0Vqpc16B+pck1hw4Bra3MQw8PM3ht2LD0BbuuNmt0Luu+ro47iooKjnoJhxmE9+0DPvpR4J//mfn1VIp5cL+fOy6LhY9DxL/r6rhT27qVAX5qijsGn4+t8qYmdrDG43wtXi//HwgsXWZtUZZeCqnpQkTweg/C5Wq5qtov8lxu9yT27GnCW2/1wmAwZumUXB5dykKdgH5ddXU1SkpKUF1djUAggEAgAEVRAAC9vb2IRqOw2+3IZDJIp9OoqqqCEAIWiwUAcO+996KpqQn33XcfDh8+DACwWq2YmJjA9PQ0TCYTotHooiN0Vqpc16BeaNhhbnx4aytbnOEwsGYNg54M69u5k0FvKSiXqw2LzEcB1dVp8fNCAK+9xg5Rr5d5bxkiuWqV5hANhXjbRIKds2YzkE7z9ZeXA6tXszU/Osq0S0UFb79tG7e9vZ3pG4eD793evdyO3HtbzEC9NrKYwlu5NV2EyMBuXwOXqwUWS+WsfYkINlvDFWeQ5p5ryxYFGzc2IxrduSBHXV9fj127dl3GqdfV1WFwcBAVFRXo6OgAESEajWLr1q1ob2/HuXPnMDExAZfLhVAoBEVRYDKZUFVVhaGhIdjtdtx+++0wGAwYHx/H+vXrsWvXLhARpqamUFVVhXA4DKvVisnJSfXeJlBSUoJQKLSoCJ2VKtc1qPt8DCAzMwxQTqcGyBI0c63l0VEGuo0bNU6aiKNDiICnn2aLdSnK7S5FWGQuBTQwwNfS0AC89BJTJWYzA7bbDfz61zzq8HoZqGVMutmsHS+T0UIbKyqA11/X6JkjR/ia9SOE9euZnkokGNj7+rhjOXBgdoJWsTzx0kuhhbck8EcinZiefgUOxxYQGRAKncDY2GuIRDphNLqWtGhXvggaoBdr1rxvwY6CiHDw4EG0tLRcFk/e2tqKYDCIU6dOIZFIoLm5GVarFbFYDFu2bEEkEoEQAhMTEzAajchkMhgeHoaiKIhEIhgcHMQtt9yC0tJSxGIxjI2NYWBgABUVFXC5XGhoaMDAwADMZjMmJibQ2NiIYDCI3bt3o76+/oZPVCoY1InICKAdwLAQ4gNEtAbA9wGUAzgB4PeFEMmlbFxFBdMD0agGKiUl/L+UXGvZ4QB++lPmmiUnLeO0h4e5g9i69eoTjqRlPDLC55TW79UWBZMdRSDAkShmM1MpRNz+eJyzS9et49FHIMDnNRrZOjeb2QLfu5e3N5mAU6fYsu/u5vvU1gb80R9pI4REAnjuOaZj9PelurpYnvhaSyIxiECgFSZTKRQlCoPBhUCgdVbooR74U6lxRCLnACiwWhuRSPTDaHTDZHLDaq1DMNgGi6UaBoP1qsvtXm0EDRGhoaEhS3cMDAygtbUVfr8f58+fR3d3NxKJBBKJBFwuF6anp/Hggw/izJkzGBkZARFl49QlN+90OuHz+TA0NITR0VFUV1fjxz/+MX72s5/hPe95T3Z0YDabsXr1amzYsCGbwLR582YIIbKJSlda4uB6l8VY6v8ZwFkApery3wP4JyHE94noCQB/CODflrh9l0k0yin9REwF5FrLFRXMoU9MsGXf28vLoRBTMdPTV59wpB8dBIPcidTWMtDu3n15xMtiKAzpPA0GAbtdu+ZEgkE7mdSoEYuFtysr405repqPu2sXA3EwyM7iqiq+L5EIO1jXr9ccuw0NwIkT3Cnk3peenmtfnvhml2RyErHYJShKOGsNGwxOJJOTWeDUW8wmUxni8RHEYn0g4uEZEWA0OgEQYrELGBl5EhZLDa6m3K4QAoqSQDI5AqPRAbO5AoDA1UbQhMNh9PfzzGyKosButyMUCsHpdCIQCGBsbAxVVVWYnp5GKpVCOp0GwIlJDocDq1atgsFggMlkgsPhgM1mw9jYGBRFwdjYGB5//PHs6KCiogLt7e144403YDAYcPToUTQ1NaGvr++qShxc71IQqBNRHYD3A/hrAH9CrCHvAfC76ibfBvCXWGJQn5oCNm3SIjL6+9nafOYZ4OhRrZ6JPkpECKYnHniAQU8f/SIE8MQTTE8EAgx6MzOzLf9CRD86aGridl66BDz4IFvI+vdnseVzKyp4/SuvaFmfkYhWsMvlYmrl9tsZdM+d45GC5NITCeDMGV5es4aPHwwyrTM1xWUDpqf5Hnq9DM5zReE0N3MHeq3KE18PshwjUL0IkUQyOQybbQ2IGGTi8fOIRDqyHLneYjaZPLDbmxAOdyCZnEA6HYDDsRUmkwep1BQSiRFUVOyExeK94oQj/cggkwliauqnsFhqYbevQ1nZ7lkRNIuhMrxeLyKRCABk48ljsRii0WiWWjl8+DBWrVqFYDAIACgtLUU4HIbFYsHv/u7vwu12o7u7G7feeiueeeYZeL1exGIxxGIx9Pf348yZM3j/+9+fpWDa29tnAfibb74Jh8NxVSUOrncp1FL/CoDPA3CpyxUAZoQQaXV5CEBtvh2J6DEAjwFY9E3zehnU3G5ellEb69dr8dQ7d+aPINGDqx40m5qAX/yCE3MAdiIeP764LMrc0YHXy8BrtV5+jMWUz5WVE/fsAT73OeAb32AAtlp5X7udt5mc5JFIJsPgbTTydqkU/zc4yJ3avfcCX/sadxQXLnAIZCTC37/6FTuQ3/te5s7z3cM9e9hHca3KE18nsqwjUCILrNYapNN+AEAqNYZ0Ogq//wii0QsoLd0Fp3Mn9DHnDsc2EBlRXn4/EolhxOO9iMcHkEyOwmKpVa3qK0840o8MbLYmpFJTiMUuobLyQZSW7p1VRXExlRYrKiqwd+9ePPXUU7BYLEgmk4jFYjCZTNmEotraWjQ0NGBwcBBCCJSWlkIIgUgkgq6uLmzbtg0f+MAHYLFYsseempoCAMzMzODYsWO4//77QUR5wyodDke2EykkNHMlyoKgTkQfADAhhDhORHfLv/NsKvLtL4R4EsCTANDS0pJ3m7lEH/I3Ps4WZ3MzAxu3ja3P+ZKIci1lv5+B8+67eQRQVrZ4nngx2aMLOVMl6OsrJz71FH8fPMiO0gsXGIxDIXYYv/02j1qI2PLu6+PfJhPTMB4Ph0K+7318/8bGOApGhi1u28YjFH2i0Vz3cKEErZUcHbNcI1C9WCyVsNs3wGh0IZUaRSgUhc22Cg7HBphMbgSDbXA6d14Wc+7x7EdFxf0AgHh8ANFoF5LJCQSDbVCUDBQliHQ6iHR6BibT4oaiuVy6xeJFJhOBwWCdZYUvVK1RD/qy0Nbu3bvx8MMP4/Dhw3C5XBBCwOFwgIiQSCSQSqUwMTEBi8UCo9GIQCCAeDwORVHQ1dWFxsZGPPjggxgaGkJZWRnOnz+fLRfgdrtx/vx5vPnmm9i7d2/esEqn05mtFaPviOYKzVyJUoil/i4ADxDR/QBsYIvmKwDcRGRSrfU6ACNL3Tg9qHR2csnZyUkGQOk0LS+fH1QGB7UQRxkBEwjwssfDID8xodUoLyRxaDHZowt1AD7f3JUTt27lkURvL68nYnrFZuNvj4e/rVbm20tKtMSjWEzr8KqrgSefZJqmu5vXK8rliUb5ErHmS9C6ASbvXpYRqF70SUKZTASKkoLNVo90mi0XIQjp9NScMedCCITDxxEMtkEIQiIxilDoOJhVAqzW1QiFjsNmayiYVy80G3Whao0S9Kenp7OFtr797W/j4YcfxiOPPAIiQl9fH8LhMPx+PyKRCM6dOweHw5GlXzKZDCwWS7YswOnTp9Ha2oq9e/di//79mJychM1mw8zMDKxWK4aHh/Hkk09idHQUBw4cuCyscvfu3Thw4MCsKo/5KKOVHCGzIKgLIf4LgP8CAKql/qdCiN8joqcBfBjMPz4M4Llr0UAJKmouAkIhBjGLhUHl3/8duHiRwczpZEt+xw6mIiwWpg9On2Ygt1p5n+lpXu7ru7xG+aZNvM184LSYEgMLdQBeL3c2QjAQx2Ic4VJSwu1vbgZeeIG59EhEs/SNRi0qpqSE9ysp0QDfbmfaZXCQt9u7V7sPQjANdbWJRit58u7lHIHmtCML2OFwJ+Lxf0EiMYl0elq1MO3IZOIIh0/CbPbC6dw+C1xyo2fs9nWIx/vgdr8LFksthFBUTnw1rNZapNNTC0bFFJqNulCCkc/nQygUwsDAwKxCW2+99Ra2bduGNWvWYGxsDOfPn0cmk0EqlYLJZEIgEIDJZEIikYDBYEA6nYbH40E6nUYqlUJPTw/27duHhx56CL29vSAinDlzBtXV1fD7/Vi/fn12xJAbVinBWR+Vk+fZXvUkIMspVxOn/gUA3yeiLwM4CeCbS9GgfMN5gOPLJyc1a3XjRrZE//3fOfJEhhO+9hpbjzINfnKSrXBpMcuqhadOcTii283AOTnJxwuF2LJtbWUr2c9UJzZvns276y3Y+SiIhTqA+nrmrh9/XIslLyvjtnm9vN+2bVxTXXLmcsRRVcVgajZrWaeKwhb5b/wGJxXpJ8rYto2dp729fN0DA1fHkb8Tk3dfQ1m2EWiuyCQhRVFABNVKt0IIM5LJPje6g/UAACAASURBVAwN/RMslioYjS64XLvgdO5ALHYGAEeQxGIXkE5PQ4gUFCUORYmCyIVo9AxisR4oSgK9vf8NJpMbdvsmEIl5o2IKzUadK8FIUhlerxfRaDTrEI3FYojH4ygpKYHFYkFzczNeeOEFlJSUZC1tk8mUpWSSySTi8TgsFgvcbjeEEEgmk2huboYQ3I82NzejtbUVyWQSfr8fTU1NqKioQCQSyY4Y5gPwfLIUk4AspywK1IUQrwJ4Vf3dA2D3UjZmruH8jh3AL3/JACvD/E6dYoCrrmaLNBxmHlkCooxKCYX4ODKRaWKCAS2pxjNIXr6/n8EylWJawu8H3nqLOxCAAf5jH7vcei+EgpAdgJxi7uTJ2eB+7718jGBQG4XIfSsr+Tud5mtPpfi3jGW3WtlKX7WKwbyujoG7poZDLXPru3/qUzx6WYop/K715N3XUpZzBJovg1QIgbGxryOR4BDCTCYOIIpUahzxuB2ZTBhWawPGx/8Xhoa+inR6BgaDCURORKMXYDI5QaRAUTJIpWYQDrcjHu8GkQ2KMoNUyq9G09TAYqlBIPAWAIKiRLMTaOgnlNBno86V8ZovwUhPU9TX12PPnj14/PHHkVItlrKyMoyMjMDr9aKlpQXbtm3Diy++mLXy0+k0eIKPimxJXplgZDKZ8NBDD2H37t2zygDbbDZYrVbcdddd8Hq9ajjmlTs/l2ISkOWU6yqjdK7hfCDAn9paBlshNGtQlpwNBNhSDQSYhunu1igYIt4uneaOwOcDbrmFO4RIhDuIYJAt3kiEQamvjzuGdeu4bVNTbC3nUgtztTm3HIE+/T8X/KenueOSpQ0cDnaQHj3KE1uUl3PbHA6NcxcCWLuWO6j+fu7MzGZus8fD1ri0nGVJ3u5u4K/+ijsCt/vqp/C71pN3L5NckxGolHwZpC5XC5JJH6anX0QmEwaQgcHgRDw+CIPBBKt1NYxGO2Kx86rlnYbZXI5MJq3SNAGk03HIsgFABtHoBRgMRggRQiaTRiYzglRqHKmUD1ZrM5LJPkxPvwCj0QVAQVnZXWhu/rtZwD5Xe/VWvp7KyMdD33vvvTh06BCCwSCsVmuWHweAStViEUJkLfRkMgkiylIwJpMJHo8HTqcT5eXl+MIXvoDh4eGsJU1EcDgcGBwcxPHjx+F2u2Gz2bB///4rdn4utnjZ9SbXFajPNZxXfSZoaGDLNR7n/5uaGPBkdEw4zOCWybBVGwxyRIjRqDkYUykGx1CIre/z52fPJ9rXp/H2JpMGdgYDW8S51EK+NucrR9DUxMfOB/6JBIOr08ngf/IkjzqIGHhXr2ZLvKqKr0tRGLT/4A8Y+P1+vi8lJVr5XRljLuc3festHsEQ8fH27mVH7OHDfK2bN3P7C60wCSx9+eLlkms9AtVLvtR7v/8w0ulJmEyVMBodSCankEwOQ4gUDAYXDAY7AEIyOYV0egZmcwVMJgcAIJkcByBAZAOQBsBAbjTaYTI5kEoFYDAIZDJp1ek6DSILEolRlJbeAat1FYRQEAgcQSj0FsrK9i3Y3kCgFX5/NcJhaxa8AeTloRsbG7Fjx45svLnD4cDQ0BCOHj2KO++8E+Xl5Vmwj0ajcLlcyGQy8Hq9SKfTKC0thdVqxdq1a1FfXw+/yofKzuf1119HT08PxsfHMTo6itraWmzcuDFvSGWhDs+FaKXrXa4rUJ9rOL9rF3PlepBcv55jrDs6+D+vl0HY4WDLd3ycAVBGeABsyVZWMthbLLzsdLJ1vmoVW6+hkDYFnMnE34AWUphI8H9SL/K1WY4YtmzR/nvzTW7bXOAvM1NlRE5zM3dQDgcfb88eBmdJJ917L5+ju5s7mxMn+HodDmD/fi3G/JVXuPOQPgZ5HadO8XdvL+//TdUeLcRRrJf5omOKcrnkS71XlAQAtnrN5nqYTF7E4wMwGEJwuXYilZpCKjWOZHJEPUYARCaYTB41ysUAIRIwGm2qvpZACAVm8yokk5Oq9W+A0WiBEAakUmMgMsFgMGXbABgQj/dcBuq57QUIL710Cd3dX4fNVp0FvJ07d+bloeUsSU6nE3V1dTh58mS2kFdXVxdWr16NVatWoaqqCgCD8Ntvv52lUBKJBJxOJyYnJ1FaWpq1ljOZDI4cOYJjx47BZDJhenoaGzZsgNlsRm1tLQ4fPpyNhskNX1zI4bkQrXS9y3UF6nMN52tquOxsZycDrd3OwHXgAEeGJBIM1idOMCBPTDDoE/GEEK+/zg7V8nK21letYspieppjuGXNmHRaK1u7eTNb/r29DLoAA/hzz82mLPK1ubn58nIEDgdz9rngf+4cdybr1gHbt7OD1mrlc0qqqaQE+OIXuX2SC6+uBv7hH7hjisd5NBIOc9tLSzUr2mrlGjAuF59/eprbJH+vWcP3Y0R1B5aWFifKuJaSGy6oKBkIkVIpiHLV4QkYDGa4XL8Bo7EMRGYkkyMoLb0T8XgvMplppFKTUJQEDAY7TCY3MpkYDAYzAAPM5lUwm6thMDhAZAaRESaTG2bzaihKBAaDDanUBGSZAUVRIEQcipJEPD4wyyma297h4Sl0dAzjttt2wGr1ZsHbarXO4qEB4NSpUzhz5gxSqRTOnj2L8vJy+P1+NDc3o6ysDA6HAzMzM9izZw9Onz4Ng8GAWCwGt9udLcM7MTGBmZkZuFwutLS0oK6uDkSENWvW4Oc//zkAIJ1Ow+VyYWZmBoqi4Ne//jX8fj98Ph8GBgawdetW7NixA0KIgh2eC0XIXM9yXYF67nC+ooKjN772NV4nQxY//GHe/oUXGMAtFgafmRmmHwDedv9+4I47uEMoL+fp4BSFI12cTq3oVUcHA+OqVWypVlUBn/88n/OXvwR+/nNOy5cZrnrAk23euZPpDoCt7Z/+dDaAO53sE2htZZB2ODiJ6MIF7qSE4I5o1Sq+hs2btX17exnQ9+3jjxDcpnPnuIOZmNBopg0beP/BQW7fli1M6chKl3IkIASPDiorZ9d9keV6V1AUy4qS2eGChFjsHBRFADAgmeyB2eyG1boapaUHUVX1CUQibfD7uUa4y7UXkcgpRKNnkUz6YLXWwOv9EIiA0dGvw2CwwWCwwG7fAIPBDbPZg1RqBJlMBEIoSKXGYTbXwGi0wGKpRSo1gVRqEqnUFEwmN8LhC9ksVsmZW631cLla4PcfhqIkMDQUhM1WDYuFE5r0QC65ZyLCkSNHcPLkSdTW1sLlcqG5uRmhUAhWqxVTU1Pw+/0QQqCkpARf/OIXkU6n0dPTg9HRUbS2tuLMmTPZsrkyZt3v9+PZZ5/NWtHbt2/HxYsX4XQ6cfbsWUxNTSEUCmF8fBxerxeVlZUIBAIYGBjAmjVrUF5evqIcnlcqyw7q+cIB5XB+YIBBPTeC4xe/mJ3qn8nw702bND5ZP4GytNgBnmBZluc9doyt7PXrtRoq69YxdSHnBt22jeu6yCgUovyA197O/LRMBnI4tNrnkp6ZmWFAj0R4fSDAVrHLxffh4kUG59pazdJXFF5W5wPIRts8+yz7AywWLYs0FOL2E/GoRnaM73kPT6YRCPD2q1fzfR4f52+XS8sDcDpXVhTLShN9uGAk0ol0eiZbRpfT8S+ipuYxOJ27MTb2DQSDb0IIBdHoeRgMNjid22GzNWa3Ky3dCyEEiAwIBI6p4D2tzvwThtf7ewiHjyEW60EmE0JJyTp4PPdkO4xQqBXBYDtcrn0wGk1z1Ish1SkqUF5uh6KEoSgKjEZjFsg9Hg8qKirQ2dkJADh69CjkXKGxWAyBQADNzc2Ynp7OVk1UFAW9vb1Ip9PYu3cvRkdH0dPTg+7ubkxOTmZj1W02G9xuN2pra7OWdmVlJWpqamCz2XDu3Dkkk0lkMhmUlJRkwx0dDkf2vofDYbjd7hXl8LxSWVZQXygcMJ8TMhxmnjgaRbZglc/HwHTnnWyZ5pYJGBxkkDt6lOmKTZvYQk4k2NoXgi3h4WEG+AcfnJ8zzwW8gQHge9+bXSLYbuf6LTIZKLe07c9/zhay16tRIokEA66cgSgaZcCWFjWgRdvU1fGIIhxmy9vp1OYgPXeOOxCPh8/V0sJRL1/5Ch/P5dJq3shJqktK+PjBIP93A0SxXDeSLyTQZmtAKuWDyeSBwcDZnzIdn8iCsbFvYmzsKZhMZWpp5VKEw6dBZIbJ5ILX+0GUlu4FACSTQ3A6dyKZnIbf/xIUJYhY7ILK1adhNldDiBTS6Qiqqh6G1/t+EBHKyvbBYLAikZiA0TibX5f1YqSjdGZmH/x+O9zuKDZuPISeni5YLAySmUwGP/nJT2AwGLLWtcvlgsPhgNPpBABMTk4iEAigtrYW09PTWVCvra2FxWLJxoZv2bIFExMT2QxTs9mMxsZGeDweRKNRhEIhHDlyBHfeeWfWmSmrOq5btw41NTU4c+ZMNua9sbERHR0dmJmZQSaTWVEOzyuVZQX1hTIS9YAq48kHBrRlIga+cJgtXtXXko0D14cRjo+zo9FoBG69VSvLGwwyZ792LX+fPs0UyT7VX1RI2F5XF48U1qzRQg57exls3//+/KVtPR4edcgJL2IxtvA/8QluqzxfKDS7nK+Mt89kuMOoqNA4e4uFOyaAOzdZ2bG9nePT16zhSTLkddx1F/CZz3CnIitVLib6pSgLy3whgfn49UxmBpFIB2ZmXobRWAaz2Yt0OoJUahxWay3KyvbB7b4rm92pVVMMIRB4Q03cuRWKEkEweBw+33MgMqrRMQlMTT2Pior75+TMc0sCJJM+vPTSZrz99q0wGAQUhXDLLYRHHvEgkahGIpHAc889h6ampixQSydlSUlJlmZJJBK45ZZbUFJSgrKyMkSjUZSUlCAYDKKyshI+ny8byiidpzK8sbm5GX6/H319fejp6QEAdHV1oaWlBZ/+9KfR1dWF73//+6irq0NpaSlSqRQ6OzsRCATg8XjwyCOPZK37leTwvFJZVlBfKCNRD6gXLrAzz+Nh+kSWpJXlZKNRrmro8eQPIywrYxqlt5cjP5JJts4Btlzr6jRqpadHA/WlCNvLZ+2vXs18/+nT2n/79/N5ieY+n5w4JBJhP8DoKI8+br2VQxS9XnYMBwLavTAYtHK8H/ygNgIIhXjdjh1aWxsbr/apFkUv+UICJb2Rj1/nUeMRRKPdAAQUJYp0ehqZTBTpdEitxjgJgDsMeexEYghENqTTY9l9FCWCTCYIg8EBk8kIi6UZ4fCpWaGLC5UEGBurQkdHAI2NQRiNhExGoLOzFvfcswG33VaLEydOwGg0zkrUWbVqFaampmA2m7MUSkNDA/7kT/4EJ0+ezIYKhkKhWbMRnTt3DtFoFESEeDyOaDSKtWvXoq+vD263Gz09Pdi6dSu2bNkCIQTa29vR0tKCRCKBSCSC1157DQCwevXqmw7I9bKsoL4QtSEBVRak2rmTnXgnTjAFMzbG4FRWxslEp08zaMnj6sMI3W62bEMh7ggkF200MsDHYrxeRq8sRjZvZpDWh1yuXq3Ffuez9vfsYZqntXV2dqfs4BYKE5RlAlwu3ufTn2YL/5VXmGsfHeVObds2jS+XoZ9SIpGiM/Ray0KzB+Xj1zOZAJLJYcTjg0inp2EwONToGA98vp8gmRyG0ehSS+wSiAwwGp0wGMwQIo14fACZTAhENhgMKZjNVWpUlxep1MSs0MWFSgKEQjWw2UJQlLehKAz6NtutCIVqAORP1HG5XPjoRz+K06dPZ7NC9+/fj6amJjQ1NS0YKkhEqKqqAhHhM5/5DLxeL7q7u/Hyyy/D6/XC7/fD4/HAYDCgq6sLbW1tuOeeexAIBBAKhTAzM4P3ve99aLxJLZRlBfVCqA0ipiVqajRA2rmTOWO/n/nxmhqmYGSqP3B5GOHMDB9rzRqmLioq+FsIBtXubgb5u+5icJVSSBmAhgYuIXD4MHcOMuRSXxemsZE7Jzlhh7S+9+3jKJzBQR5BLDQS0E8cIudtlZN9tLUx7SIdyh0d3Gnt388dzNGjWgeayfA+cjLqIt1ybWQhekOm4+v5dSIP7PY1SCZHIUQJTCY3Sko2IpOJAaDs1HXhMMd7C6HAaHQDsCCTCUNRIlCUBIzGUgA2lZePQVFiYFAu3GqprCTY7RtRWloJIcIgciIcLkdlJamZzfWoqNiFzs42lJXx9ekrIU5O+pBMVsFiqcHgIKmBEJeHCk5NTWHTpk3ZJCWn04lgMIja2lps374dHR0duHDhAiYmJgDwMcrLywHw6MBoNKK8vBwejweBQABHjx4FEd10VjqwzKBeKLWRa9HLuuCVlcBtt2n/KYrm8JNhhFu3MuUyMcFW+m23aVEhfj+vb2piTr2xkQth6c8veX8Z5+1yXR7Drb+OyUlttiUZvaMvqiVr2egdubLTIGLaZM2a2dE7ufdCThwii5MFArxOVm/cvp2v6eJF4P77+QNoHSgRUzgAT2QtZ5FaQSVzV4wUWvEwF/wdjm1qKd4EnM7bkclEEA4fB8BT17F17obFUo54vB/pdAhChOHxvA8GgxmRSCcMBpl1OgxFiUJRMnC73w2nU0uYXagMABtehLa2ChBVZPVTUXj2rPZ2AtFBCNGC8nIfPvIRLxoa6tV9G9De3lCQbsuaLW63G+Xl5VAUBYFAAF6vF4ODg+jt7cXWrVuzU+F1dHTg4x//ODZv3oyjR49mQylPnDiBjo4OAMy7r6Tqikslyx7SWEhGYj6Lfv9+3re1lS3WcJiphmCQaQVFYQfjgQNcR6WzkyNObruNLVhZcjcYZAs6HmfO/cKF2QA3Ocn/ycqIMhlINRhmdUb19dxGGdqYSnFH8pu/qTktczsE2Wk4nezgnZzk754ezpjNBdq5Rje5lrjbzTHvW7Zo++tr08/MzHamFpONro0UWvEwH/h7vQcgBDAzcxip1BSSyQk4nS0wmTwq+AtUVn5EdTAeARHgcGwBABiNdoTDHbDbN6hO2Ajs9lthNJZhevq5LGhLzj8abUQkYoDDoQBog9PZgsnJhuxctjt2AD/6EevN9DTwP/4Hj/LuuQcwGgmK0oD+/oZsAAOwON2eLzX/5MmTMBqN2LFjB9asWYNwOIyZmRns3LkTDQ0N2f1CoRA6Ojpm8e4rqbriUsmyg3ohks+ir6tjS/j4ceaSq6q05KSdOzWqRlIa993HQNvWxhaunMPzjjuAn/xkdrihHuCSSXbQrlkzOxno5Ze5Q8i1wL/3PaZI0mnuKMJh7lxklEpujPvEBLdxYoL/dzrZyp9rRqa5RjdAYVRWQwPv5/HwPQAWXzJ3Jc92tByir3g43za54G8212J09BtIpycBEMzmVUinfYjF+kEk4HK1AIDKzzcjFOoE14HRprxzuVoQDJbA4dgCLvA1Ow49mfThwgUDenoMqtFiQHOzAZ2dPrz1VsOsoAOfj987aeCMjrL1LZPVrka350vNl7y9EALl5eVwu93IZDKorKyctd+RI0cAAFu2bIEsNnYzJBvlyooAdeDy+uWHDrFT8NQp5pcrKpg+6e3luijHj18OcAcOXA6GJ09eHm6oV06L5fJkILebnZF792r89Cuv8Ajg1Cl2uBoM3IlMTTEAVlRc7ggWQtvPauURhtnMvDzHJucH2rlGN4VG6VxNydz5fAxAEeyvRnLL3Y6MPInx8adgNLpBBDWCxgiP510oKdmMUOg4RkaegLTsAQWxWB+4HgxPeWe1NiIavZSNhc911Pp8XoyOKigv1xydfX0Kuru9WLdO0++XX2bjRq8zcrIZSQNerW7PlZq/UIEtud9dd92Frq4uNRmLsKjqinNZKivQglkxoK4XOaxzu/nj8bCTVFrbXV3zx7/rdWYhgKus5CxNOR1eSQnTMUIwrSMzR6VTcnKSI188Ht52epqpHpkIpLeeBwe5c5BJRAB3AuXlTNuUli4uqzMf2OfTyaspmTtfqeF8HelCPP0KfGfeEWFa5E2YTGVqlItAIjEAu30tLJZqEBFCodmhkvF4P7zeB2AwWLM0TyIxiPkctX5/PYLBXSgr02ifsbFdGBysh82m6fe5c/yczGYOBHA6WWeHhvLr0FLq9lxWPAAMDAxk/6urq7uy6opzWSoHDsxdL/s6VuoVCeoyvl2f3k6kVTAE5o5/lxNV6Gmc+QCuvp7pE30ykM2mDStjMbZWqqr4c/48O2CJOOXf6+VaNbJ2jP75+nz80qxaxe3o7NQqNr75JvsN6uqu/D7NZ1Vfaez95CTfg6EhfrFlLPxCHeli23ezA3sq5cs6OqXlKQSQyURhNnvnDJUksujWAxZL3byO2spKwvj4QTidLTCbfUgmvXjrrXpMTBCCQU2/zWbtE48zsLtcXIJi27Zrr9u5VvxcU84dOHBg8dUVBwbYESZrdkh+qLp6RSr1igR1aV2XlWkJRjKsL1/4nrS+Kyrm7pDnArhcDjuR4M572za2YJJJBvfbbuMp9jo6WHntdi1CZf/+/Ek9Xi/rUGMjdwYAX9OWLcC73qWB55XSgfNl7C42U1rG/X/jG+yclslaDQ18jUDhHam8vyt5jtNrLWazF0ajE1ZrIxIJjvjIZAIoLT2oi5y53AIPhdoRj/dBD+AVFQfmdNRq0S0NMBgaMDPDutvSwlgn9XvDBgb3sjJ+12pr+V1YtWp28pqUa6/bc0w5t3MnGohQ0GGF4Jf4r/4KeOMNHmLbbNzo8nL26M6XHZnPIr8OlHpFgrqePpBOGn0RLyC/9Q0UTsvMJdPTTLPs2MHO0+FhVu6GBu40pFW/fr1WTGtoiDuCsjLg7ru1YmH666iq4mOvXcuJVJLnv5LkIKlrR47wyyM7KKmTk5McZlmIMSH1/i/+gpO+pBM4EuEXvaMD+PjHF9+RHjy44uc4vaZitdajrGw3AoFWEJmhKBF4vQdRXf0oZPXEXAvcZluDWKwHdnsTcrNX53LUErFRU13NGJZMMr3Y1MTvlNTvW25ho0nOJ6AorLPxOPDVry6HbvsQCmkjFYPBAAMRfE8/jQZ9FmA+xRaCe6wf/hB4/nnmlhIJflmam9lKNBq1mWbycbNzWeQSXJZRqVckqBcS355v/cmThd1vfQcsy//KWHO/nzNZ5dDR6WTeXSY3lZdzdufOnbzvv/4rx/OGw7y+qYnL+v7Wb80u2/vLX/I5p6a4nTJ0UtZkKVT0uhYKMc+dSnEnJF/GZLIwY0Ie69Ahbp/FolWRHB9nw2b9eo6LHx7mofobb/B6YOGOdCXPcXqtZaFQyHzrk8lJxOP9l1Ey0ikqpRD9liM5qd8AGzF9fdwBbNjAxsrnPqfpttfLIY4PPcRBBNdWt704flxBKqVgxw7uwJRAAN6ZGU4+mUux9VEWb73FITxyEgK/ny+wspKtMUXJr9TzWeTV1cuu1CsS1IGF49vzrZ8PRKSiSyu2t5c7a5l1yfG4DOajo1wvXSb9uFxMtdxxB1us9fWs8K+8wsZANMpcZCbDRsH//J/87GVt9Bde4LBKaSB4vTysvRLR6xoRA3pHB5/f5WKdtFi0Do3joHmI3dk5u3OUxzKbeVTqcjF4p1IM7rEYb/uVr3BGrsHAlhsR8KUvcRtOnMjPwft83BnM5c8oOlAXDoXMv15fIExBMKjA7/dmgdnn0/TbYODnMjTEFKH0Mc2l3/ffr9XjP3cO+NnPeMRmNDIl4/OxHr3+Oseh//Vfc4jvz37GowCTifFyaXS7HqnULnR0tMFsNsDlUrBrzRrU62enIeKGHTrE0Q6ybkdrK69Lpzk0JxRia4yIL8bjYafZM89oSg1wVTybja38ri6+2LKy2UWWLJb8Si1jsGVPClyz6nkrFtSvRPJFfaxZw8/mxRdZ0aV1u3WrZt2++SZPaFFby8+vpob32b6dn9PkJPDqqwyO8bgWBdLezuusVq0AWTLJ6//0T3nYm0zyJCDhMP+W0+vdfjvr19TU4ops5VIaO3YwKO/bxyUQpJGhKNr8pX19/AKnUlx3XdJYsiKkzcbX5XDwbznJ96lT2sxKzc3sUxCCOdSTJ7ntL73E11tWxu1pbOR76PXOH3Ovz7Kdmbm8Rn5RLpdcSubCBQVnz+6Cz1efzSCuqeHnIadavHSJ6TVF4fIT27fPr99OJ+ul1G2Z65BM8rNPpfg5vvwyG8InT/K6TIb1x27XkvGuTrcJO3YchNncgn37fLjrLi/qhQA98QRfDBFf6Kuv8gtstWqT8168qE0tFonwwYNBrcaIzcbWSHU1D0Xtdj4GwOB87hxvGwiw8uuLLFVW8kuXm1Qjo2j06dyLnTuyQLmpQD03nb+9nS2Ijg4NyL1etkpkiOTAAA9Hjx9nS1aCUVUVj9hkZ2s08n6HD/P6LVvY4pEWraJArZXBBsLgIPCd77CVJKehS6X4t9nM34sNaQTyl1RwuRjQ5ahFdm4yjtjl4nONj7Nx0tPDOlpby7qtKPxODAxow2ynk9s5Ps7vxfnzvG7DBn5XvvUtPt/x49okIIkET0zyO7+j0Vf5RlQDA6z/DQ1ap3PqFHe6+/cXo2PmEj0lMzjoQ2urF1VV9SgtpSxeEWl5FlKPx8f5+XR0aKG7+fSbiP00H/wg64ac0N1s1kK6UynuEIaH2UBKp7X1Mhlv6XSb4HI14K67Glh/JEBK7rG9nXuSTZv4AD4f8KtfsSJt3Kgltkju1Grl/wYGWKl7etgZUFfHyyYTK/PUFPdU0lnwxhu8/sABzerWK7VU6MZGPpd8GNdo7sibCtQB7X4DDBZNTQysZWV8791ubduREQb/ykrulGMx7uTXrWPLQ84FCvByJsOWj6KwHqxfz8eVpYIBzYKWk2dIusJq1UZ64TC/GPv3Lz5KpdAiaXL+UjmN36lTWtbrxATw1FPc/t5e5j8lWCsKW1ovvsj7VlZqZRnGx/n+hcP8HsViWhKWx8PvVCrFHcWhQ6zH+hEoMNvBK6cnlO0qzp26sEhKJhRqQDLJ9y0c1vROSirFemC3a1nSisLP9Y47uNOWTnY51JNzDAAAIABJREFUyhwY4GcyPMzP3WbjZ5xKse6bTKwrgQAvE2nnNZkYB6PRa6jbeqvtyBF+eYNB/j8SYWWSnt7OTn45a2tZwfbu5Yu5cEHjF41G7pXGxvi7spI7isFBPt/0NG+nKPyST05qkx/rlbqzk1+OsjKt3jVwzeaOXBDUiagewHcArAanrT0phPgqEZUD+AGAJgB9AB4SQviXpFXvgOiHcurkLAD4GTU0sNUii4Bt3szPsaSEgWzPHn72Q0P8rPx+tkZefJGfP6k1Y5qauDP3+bSKkJLPTqdnW9PyJZOW1G//Nnfsg4OLo9wKcSJLkaGI+ggZv58NFCL+3dTE74GMw4/F+B5Fo/y/0cjLkQjrclcXbxcIaCWAUymt2Fosxnr/5S/zaNhoZHDYv587l/Z2jQIbHeU2RSJ8vHR67izbK5EbVbeB2VatnKIQYGpFOsqJGKDdbr6fk5M8exgR012nTrE+d3fzcwNYJ557jg0ap5NB3mDQsC0S0UZz+sQkaelXVi6xbtcJkN75UlvLF+XzaXHOU1OsnKOjfBCPh4fSFy/yC2k0cmMvXdJ4UEVhBZd1uSsqWAF9Pu4JJd8UCPC2djvwN3/DXGpt7WylPnyYKZeREa08bDzO+/p8S+5ILcRSTwP4rBDiBBG5ABwnopcBPALgsBDi74joiwC+COALS9ayayx6pfd4WMk6OrS0549/nJ/Nz3/OHez0NG83NcX7l5byyK6sjF+GyUmNl5RD2UuXNAtFWi2pFOuI08kjOYC3Tyb5t8HA+vHd73L4l+T9W1pml+ydT3JHfzKCS0/xyQCAvj7W50SCr2d0VKMZx8a0rF0JztEoUygyMUXeQyFYT8vL+R6+9ZbW/vPntZGOw6E5QGXnYbHw9vX1PFNUYyOfq62N30khNN+Vw7H4qIl55IbUbWC2VUukVS+VPsFPfYqfYWsr3+tEgulHGcVnt/MxTp3i55VMss5XVTE2BYP8LOR8BFK/5ehg9Wo2eGTp63Saca6kZAl1OzesUEYijI1x7zMxoVlMktg3mViZw2FujBB8sefOaZMKyzAx6XgSgk/8ox/x8LakhC8oFuPzOBx8zokJPnZ3N1+sVOr77+fj9Pez0sfjfGNSKT72XXddXZZhjiwI6kKIUQCj6u8QEZ0FUAvgQQB3q5t9G8CrWEGKL5VeVnlMJIAHHuAol1WrtJFTIsEd7cwMPxeAnT8+Hz+PlhaedejIEX6+FRX8vKRVLp+b1BODQYsmaW7m9dPT/JxLSlhZKytZX1wuBt3XXuMRnIxemSuefK7SFbnhtI2N/J+cjzQa5Y/fr1lkJhO3bXKSt5HZu6GQRgkC2jsg25BO8zZydOJ2s7V/+jR3ZE4n670cnpeW8vaRCB/XYtEcdrEYO6gdDu39W0q5UXUbmB0u29XFlEp5OWOSBFCA/X9PPslUW3k56/bp0/yMLl5ky76hgfVfUVgvJA7KiC55PrluZoafXV2dFuq6ahUDuRCL0+15o6BywwovXuQLamzkHZub2WKRZVXtdr4B6TQrm8nEyz09Wi0DeUF6SSa1GeXl0L2ykr3B0uEkrQ8hNB41V6nf/W7uJTMZTcHlJMRXk4mVI4vi1ImoCcB2AG8BqFJfCgghRonoCgOVlkdk0oXPx3ogJ5sYGOAXQSrOwYPszP7611nZR0a0kMOJCaYJdu1iSzIU0uJ80+nZFrg8p0m94zU1fJ4tWziO99IlfgkcDgZ5aRH19/NoQFIzr7zC++mjBubLTM7V+0yGHbQ9PQyo09MMvNIwiUZZ9+VxUym+NrNZS0TRixzWy9+SMgH4HKtWabxtIMDnkx1dMsn/SWebxcKdXV8f67fJxFRXff3sCUEWGzVRiCxWt4noMQCPAbiuKwDmq8ejr+e/dy/r2eHDrN8XLrB+j44yTo2PM+4kk3zvZZ5DIMDPS9/RyukhTSaNmpP7Wiz8DOPxwnV7wYx7PYcqBDtrwmFW1HhcqxwmQxLlUDmd5k8gwC9trlLnk7Ex/ibi4bvTyaCcTPLNksNVgI8J8EVLC04qtbRQZIwpsORTkBUM6kTkBPBjAH8shAgWWnT+elb+oSF2EO3ePX8SjowAsdsZ9Hw+1p3xca3mTH29xgUnk/yRlosUIgZVq1V7zmfPsg5KR6vZrMV0A7z/5KRmyc/MAE8/DXz2s5fHk+dL8MkNcQwEuOMxGlnnLBbN0SlpwkyG10ujQ84FW4ju6yUa1ZzEknOVI2GA3wk5MhCCLcV4nLe5eJHDMHt6Lp8QZKnzOK5Et4UQTwJ4EgBaWlqWeAyxNFJoxroEVTnX78AA66ekTvx+fvbSpyHzHzKZ2bWWJEsho7vKylh3ZAeuKFoJXmBh3V6w/XoOdWaGX0oZoZJOa5yotGbica2hMgNKco2FirR0gkE+l8PBx9Vb+PKmpNNsOeVT6muYnFQQqBORGaz03xNCPKP+PU5E1aolUw1gIt++17PyF5KmLhWroYF/RyIMujJyQOrKhQuaw1w6Q6WBoD82ESu70cgx3JJjJ2IqbnKSjzk9zd+Tk6xHNTUa7dHRwS+etGjmuw6p95OTbGyMj2u6GAxq+ibjjQGNI81k+Fr0UQyLFXkM2VHMJek0+zDe/W6+D5/8JFvp+YrkLTZqYj65Gt2+3mU+vZD5Cp2dbCFv3qwZrnKeAJttdhiuTO93u7X5AiTFpx+xEbEOjY1pWDc2xscnYt1Op9laTyZnUyp63V7w/ayrY1L+2DFW8GhUC1WT9Id0fgG8rI9YWLeOHT6SeilU5MUqivaCzLVdTQ1zut3dwKOP8vKPfsQX6nZrYZhLqNSFRL8QgG8COCuEeFy36icAHgbwd+r3c0vWqndICskw1YfX+XyzHZ/ptDY5x8SEZt2WlfE2Mt7cbNasXquVa2AMD2ujQUnJyHjehgamZe6/n7d57TXWgWiU9eDixdkWzXzXUVfH7fzOd7gtshOS+iRfzupq3reigoE/neZ2yQm6Zbz9YnntTEZLqiotZcCQlroEBPlOyJj9D3yAAX1oaO65XZdCbmTdBubWC309nvFx9hHKQAyHg78B1s3ycs3AADSKpaxMM1LlfLmS6fB6tQ5CshKyQ5eOblkQsb9fy+uprZ2t2/OWkRCCe/zubr6I7m5+SWWmn4wT9ng0Z2R7uxbfmUzyyfNx6IWIdDDIGGg9zwpoSjo9rSm1TNWVQ6Py8muSUVeIpf4uAL8PoJOITqn//RlY4X9IRH8IYADAR5asVe+QzBX3KqND9PVTRkdZYWU9aaORf0sHoxyBSUeQzMR0Orlzrqnhc+zcyceanuYhrX4oK8McpZFx9qxmFPj93K61a3nfnh7udBoa5o/fHRxknVq9WnPanz7N3yUlGqVit3M71qzha5ajDjlakNeWK9J6ki97rkj9lSGb6fRsMJcdodnMbXzggbktdD0XvERyw+o2MLdeABqtUVbGFnNvr/Yc3G7WQZOJcVJSiQDrQSzGz622ljvbmhp+dskk604mw99Sf9JpLexx9Wped+kSRzo9/7w2emxo4PNK3Z43Ll0OocvKeMXGjRxtImPSHQ4+eSrFPUxnJ/+227UXT/Ze+UT2QvOBvuSmUqnZ/8shucPBL+0DD/BNeuKJ2VxSf7/2giyhFBL9chTAXGfdv6SteYdlrpjufPVT2tq0RIvaWm1+U2lpm81a4oVUYjl1l9fLCvyhD/F5/vZv+WVxOLT6GjK00uPhRJ/+fl5/xx08CjhxgsHVbmdHututDUPni033+fhlczg0q8vp5OV0mgMELl3SEkM6O/m3zASMx2eHW+qH2XLZYtGsqFyRL7PNpnH38t7LbG7Jp2/cyMcaGnpnqpfeyLoNzK0X+sJ2Hg9HJ7W1sT5nMvwc4nF+DyT+ORwayCuK5tBPp7Vp6u65h9d3dvL55QhVGkIOB59f1p3Rh5Unk3y+Xbtm6/acOReSm4lGtSGvDMfJZGbz6HLqJekklWAvxWS63GqRDi69A1QvRiNflARlq5X/l1EAMuGkvJyPIytHzsf1LpHcdBmluZIvTT1f/RSTifVDznkquejhYX6eMjLAYtFicgEtnvvb32Zece1aHmJKB7zUG8nFO538IpWV8TpZ5ndwkI/T1MTxxIODs30r+a4D4G2sVk03zWYtjBLQLC+zWXPM+/1aeV3pyJTnyAV2QEssko41KQYDD7HlOxMO832RmYczM9pxy8p4GL5tG7/Y75D+3/CSTy9yaY1t27R7fOwY64HDwdu1t2sUmuTQZX6G1Pe6Oo2uFIItb0Dj3zMZ1mnpM6ypYX3r69OixKRu5voN59Lr7EVIp6iMTZdgLiMbSkr+T3tnHt5mdeX/z9XifY1tQojjLDhOQkIgG4RSugCZUlqgAyVsLTBsLTttaYGW6VB+TIdMWTtQQqAUWgpMGDoQUpZCp0yBCVkI2SAkTuLES+zYseNNsrXe3x9XV5JlyZbkVfb9PI8eW6+kV1evjr7vec8951z1o9I9JyKNGnpeZmpvLHyyKdxb1144hAqU9CVnWZkaR/hST2vWqA933319xJIGl3Ev6tGI1j8lL0+1GX37bSU++urO6VQCqWPo2qPJz1ff7eTJyoamT1dpYw0NKl6uY5BOp5qkOnJEVfMVFYUKdDo7lVdTXR1a73T7dnWiuOyy+OZWpkxRhW26MlD3atFhEK9XORMZGeq2aJGK3+/dG8oK00Tadnp6z0ZOEEr91dkO+nGrNRSG1FW8OswUnvHz3/+tjs8w2f+4JFpY48wz4bzz1CIo2r7b2kLOZmOjshWPR80j6ewYjyeUENDcrF53wQXqakBXZOuc9hkzQh77nj3qt3TokHpMx9m3b4crr4zDtvWHWL8+VMkGoRhpa6sS546OkKfldPYM9EdDe/2a8ET88Bnh9PSeHnp2ds9cdZ3vKWVoJjrZNSQTxIh6FGLF8qZOheuug69/XdlRd3cofVVXE+vJRgjlfefkhLJc2trUBKgQoUUJjj1W5XNrz7SjQzkgDQ0hwy8pCS1C0NamLknjCcUJoXq3L1oEDzyg7GvqVLVv3ZclvBK6vl69X1dXqEAoFrpIT4dRwkMw2qHxeJTd62I+i0Wd8PT6lzo/Xdd3bN2qWrUOk/2PS/oK14Xbt64udjiUsDscobYAhw4pO9ROTHe3+v4WLlS2umyZclQqK5X9zZsXSot97TXVP+izz9R8VVqamkfRWhyXbesPMWmSmnzSWQm1tcroLJaQ13TUUaFeFzpbQIt0JLEmh8Kfq3+oOhdT52wWFoYu061WdXB0eGbdOrj00uTWkEwQI+pR6K9/ir4c3LxZVcmddFKoK52uwpNSZc7oiUhdmNTWFpqMP+449fzrrgtle+gl86IZ/rx5SngPHEisAEeHTZxOFb7Rczvaw9JXm9r2s7J6inO0uaJwp0U7JBr9u7BaQ7229f+6d4y+3NYnlezs0JXsnj2qjH0Y7H/cEiusobdPmRIKA5aXh8J277+vrig//TTUYqCtTTkuuhZBX2Hl54ccWSFUEeUnn4RsYfJk9V1rWygtVbYWXpcTlfAy05YWNeD2diXYhYXqEkHP9BcXhy4NdQVn+GLGiaJjifpD6VzdrKzQvvUklI6v6p7rzc3qrBdrKbxBMnAj6jGIGcsLoFvJ1tcrezn+ePWd7d6tvrcJE9TE5iefqO8TlAhnZKhLzMJCZVff/KY6CehYfV+G73CEXhdPKCLcburrQw3JdHWm7iYqhBL09PRQTxrdpS9S0LWD019qo86W0GmRLlcoHVRfurtcPQufrFaVHVFSMmz2b4hCeCVne7tKeZw8WdnF5MkhkW9pUd+Lrnvo6FC2pZM6dNvwDz9USx0uWaKuGLXwFxaGmufpNhxRr8jCDaCoqGeZrC7hLipScUo9M5+Xp0RcN/NyOkMVoAMxHt3aVHtrekJM/zB0GpnO3ZUylJerf7T9lsoODCPqSRBp9K+/3tPoi4rUd3PqqeoyVXsu6enKqyksVJekxcXKPleu7PndJmP4kcIX3pc/lu0LoX6cEKqStViUoxM56Rn+PvEQuXiCFm5QJ7uMDHUiOXJE3dxuNeZTT1XHcZjs3xCGtqHt29X8z7x56squuVk5Fueco0LYOhuvulrppZ6It9mUvX//+6Gr0nnzQvMpGzcq2168WO1ft1O58so+mnpFGoA25PClyBoalMeic9JtNmXQupxV5+XqXsGRxArFREM3doJQiMduVzftKR1/vBpTc7PaVlCgJtUmT1YHTVd8HX98z4MzSOldRtQTREo1kbR2rfJYli1Thh1u9Po7379f5QDrLAGnU4lqS0voEjNa6l6ihh9N+KZNU++v9x1p+7o3h85M0NlfEyaEMnOSRcfSw7MadLaZvp+To4Rcr3Xp96umUm73sNq/IUC4DTU2qjC17jtVXKy0qqxMHf9XX1Wet14Ew2JRDsfs2aEw46RJSsv0yT08g0mHmXXmn14sKOpJOrJXgJTqsla3U7Va1Rt/4Qvq0mDfPvUjrKpSIq4zZHQ2S6wPn8iB0mjPS8fRdSvSujp1ALTR6xajV10Vqtj6/HP12gULBj29y4h6AmjDf/11FWapr1e2tmBByOilVLHGgwfVWo16ItHtVoaflxf6/iB26l4ihh+tR4bOYND7jrT9v/0t1FLVbg+9t5701K8Lj5sncpzCx6mLq3QWjMejfvg6tq7XwMzKUmMeRvs3BAi3obw85aToSvapU0MhvwMHQt089YpWubmh1s36aixWNaju5a6X0+v3JB2ZX5yXp/7qHsI6k2X+fFUMUlMDv/51aMWb9vZQkYQOmSTimWvs9lD+ZSS6NHvGDPW+Ho86gFlZKvXMalWTY3v2KE9wypRQm8pp00Id9QYpvcuIegJowy8vVw6BlEpwpkxR30lJiRIevbivrnvQedwuVyierL+/wTD8aD0ysrNDIRS9j3DbLypSJyY9l6OLSkCNU+9PVwRCSOz7KrILPxGFh2/0/7pqNjdXjU+36j36aPjqV5WQ7Nw5bPZvCKBtSAgl3O3tymN/4w31PVx2mXrexo3qe9i5s+e8i9+vXjNtWmjx9cgMpsWL1VVqtPWaY56kw88OemIyIyO00G5bmxJTLbZlZSq2uWZNKL9WV1XpCaTwlEYhlIHpRPxYRMtZj2TCBHXpWViohKGwUG17/301jvR09cNsblZisW+fmkDLzlbtMgepp7oR9QTQHnRNjbIT7elaLHDjjaFq1Koq5RVXViqB1OENv185DEuXhmLig2H4kV6Rfq9Jk1T4Ii+vt+2ffjqsXq36Gdntyt5cLvUj1c3sdDwc1P0JE9T9lpbYjo72xvtyhNLSVKglO1uNqbFR/V9drUR9GO3fEEDbUHOzOoHOmaPsb9489V3qpQctFhUay8tT9pGVFSqenDtXzXVEq3LW85tvvBHqTDptWmi95pgn6fCzw+7d6hL4mGOUIRw8GCptXbkyNNmydKlKuv/rX0OTmLo1pJ7A0ScJbfz6h9xXVkysx3RWxZe/rD6onkTat08J+qFDahwORyh/d8qU0KpHWVnqua++OiiTRUbUE6C4OLRu5nHHqe/o0CElMDq3Vns8kycrEdT57DrF6/LLVQdC/b0NhuFHrnKjsw501V802//Wt1TTuKeeCvWI120BMjOVuGdkhIrjMjLUvnUufl8OS1+PWa0qL/9rXwvZ/7ZtyqZ37Rp2+zcE0Db0+uuh5Q1nz1ZXi9XVobVk/X4l4CefrGzUalWe/CmnwPLlPXtThWeQ6bWX581T9rF/v/reLRZV+BSzBiE8H33VKjXhVFSkBvT662oiSw8s/HJ2xQoVg1y9WnladXWhtEct6FlZoUY1/aGzaqI17tIru1dWqh9TTo66JG5oCHXB0yXjen3Tujr1/qefPuiTRUbUE2DKlNDajXrO5oQTlHjrvPFoht/YqF53zjkqJ12HSWBwDD/cK9q+XZ149MTi4cOxbf/880OvWbdOOUJ6IYyOjtCaorm5oUK7/kQ0vJVAZC95nWl2/PGhjow6JDRC9m8IEK6dTz2lhDrc6dT/FxWplZQKC5WdTp8eX6NB7exYrSpEOW2asrdvfCNUjNfn4NLTlYeuPRvd80X3o468nNXZAZmZ6sezdavyBhoaQhMBujlNbm6osipadgyE4qKRE0+FhWpMDQ2hxXqPP16dfKqr1b4//jg08+xwqM8xebK6fK6rUycBfUk+CJNFRtQTQAhlwFVVodah+fnqu9PtekfK8PXJ4fBh9b4666A/29dNyz78UL1nZ6f6LM3NynZ17rpOr3S5Qtli0dAZZOHoPty5uaGiOz3O665T7zVC9m8IQwgVuaivVydNXUOwZImaI9q0KTRnUlSkUhfj7RobGSIsKFAFnscfH+fVVuQO9Io11dXKOCZOVKGV+nq1TV/ehf+gCgrUY/qyt6lJeVw6gyWeVWB0L2q7XXk+drvaT2urMuAJE1Q2hY5bFRYqgW9qUmLxpS+pg/b66+qy9KjAolpTp4ZOEAPEiHqClJWpCW3dlretrXe73hEz/Cj70NWheuHhaH1UtO0XFSnb0i1/daM6PS/g9YbmnHSFYSz0HFVxsXqNXtkrLU2FWTIyVJgoJwcuuWTE7N8QQbRqaikH3jW2zza6ie5ACJVN0tmp4uagBjllijKcaJVOehZY59TW1YWaNuk+B33F08MnivLzQ82TdPhGt7Ds6FAnjvBeF42N6scyf77y8HQS//TpPVe+iavpTf8YUU+QWC0EoqUVDrvhR9mHx6NCPx99pE4SOTmhCVOdpaJPBFIqh8ZqVVcjRUWhxRO03UMogyUW4fUY6emhBbhbW9XvSDfCy82FF19UAj1C9m+IghChY6urkcPTXJNJK+2v9UZCO9i+XcXvpk8Ptfvcu1cZW36+MvSNG5VHZLWqRYBLSpR3IaWa1NHLLum1HfsyaJ1O5narfevJJQj1mfF6Q43DMjN79rrYvl1NlM2bp17b3q5eqyfiOjsTaHrTP0bUkyBaC4F4lsaLZ78DMvyIfTQ2qq57e/cq29u6VQloXl5owjTcM16/PrTMXUZGKJzj8ymb0157X4VJurd8SUkobDNhgvrtQChurluxTpig3veee0bE/g1RiCxma21Vwj5lSkj/kkkr7a/1Rtw7OHxYGZnNFqpi02sm6v4vu3fDP/9zqDVAU5Oa5CosVLFRrze0pmlGhvoB6F7AkejwjC6c0J3wdKGF9nj0pJMuida9LnQjnfADevTRoeos3Z+436Y38WFEPUkiy/KLiqLnnA+74Yfto65OCfmMGcrmurpCIllWBs8/r6o1Dx1SjosO95SWqt/E9u3KlnVHUb24hl4oOxp+v/q9dXUpO5YyVNiXl9czFp+WFurWOEL2b4hC5FVnaakS9fCeRSPaNTNykQC9SosWeZ2XqS9LQW2rrFQCf/zxocUQdLFGbq7y4JubldhnZYUWA+juDnlXulGYzqZxOtX76+XOtMHX14fatcbK7xyiNqRG1JMgVj+SxYvVZNJoaRe7b1/oikE7GEIop6WgQE3YT5qk7HL6dBUaqahQse5w2z9wQO3L4wnlI9fVKdsOX5IvPT3UG7u9XQnv4cPKCdJzWU5nqH1GWlqof3z4ogjDaP+GKERedepq5FNPDfWlitaqYtiaroUvEqDX4TvqqJDhNTSEWkTqQehY5IQJoaZekyaF+mSUlCgxnj1bhXZKSpRHkZ+vMhYmTFB9sNPS1Ov8/lBBkm5FmZ6unnfMMb0/fKS3VlY2ZG1IjagnQbT4+caNalJ0yZLR001wxozQFYP2lqVU9qpTc7XQ6x+x06mckXDb188vKFD7zMoKve6oo9RvoKNDHY8pU9QJoKtL2W1OjjpeGRnqt6Z7Lun5qn37VHZLeDHRMNq/IQrRSvylVCf5aFeQw950TYhQPu6OHWrb3Llqu27a/oc/hMq+QQ0qP18ZWmmpMnS9/mRBgfLYGxrUTaeEud2qoEjHD9PT1YnjyBEVNy8shLPOUo///e/KkHV8s7297/7Yg3FJHgMj6kkQK34e3i42GsPdQvbkk1VNxPvvhybqS0rU/zqscfTRak5H/4izspQAz5oVsv3MTOUtO50qTg/KeVmyJLSC16efqtd0dKhQDqjflv6Mublq/V2nU+1Lr4OpW1DX1sY+bkNo/4YoJDphH8vJGdI6AiHUG0aK5tSp6odWX6/ii1VVavvRR6sqU4sl9ME6OtSP5FvfUhM7q1apH0x4cZPOJsjPV/uYODHk/fzjP6rn/s//hBYD1lk4bW3RY6/DIAJG1JMgVrOivuLnI9FC1mKB++9X9rpvnwqxHHOMCrMUFYVyj7OylO3rBbXPOEONSY+1sxO+8x11wvr0U7XvmhoVtty5M5Shkp+vHtNJCF1doayxc85RFaxvvqmyV2y2kD07HNEnlE0f9ZEh0Qn7wUgSGFT0B1i0KOTJz5sXGkz4BystVR7Fvn0qDqj7ZhcVhRYz0An73/mO2qcur9W5nnr9RV0taLWqH1G0/tjDIAJG1JMgmdTDEfFmUO91yinqppk2Tf0tK1PjbmoKJRDo1r4Q3fYnTVIe+AcfqP3oibSGBvX7yc9XJ4np09Vzu7qUzV59tSrtf/ddFbM/eLDvFgimj/rIksjVUTJOzpATy5OH0AcLN7KODnUJ6fEo70VKVVZ73nmhggt9ZtP73Ly5d7VgZaWqFIxWLThMImBEPQmSST0cdd4M/f9wo9m+xaIEWXfkg94TaXrSXzs04Xn88To1I3USNCTOYNRXjAjhRqZn+7dtU5M9ublqncqlS2P/sKNVCx51VOxqwWESASPqSZJonHdUejNxEimw2dkq3Bje5CnaRFq4kxStBUJfTs1InwSFEGcBjwJW4Gkp5f1D/66pyWDUV4wIkUa2cKES9FNOUbH1/j5EomezYRIBI+rDRMp6M/S2/Wjhxv4+S6JOzUieBIUQVuBxYBlQC2wUQqyRUn429O+emqTkZHakkemuqmGiAAAgAElEQVTmXqedFt8HSfRsNkwiMCBRN95M/KSsN0N0248VboxFovY8wifBk4A9Usp9AEKIl4DzACPqY4nBMLJEzmbDJAJJi7rxZhInJb0Zott+f+HGSBK15xE+CU4GasLu1wIn9xyfuA64DqAs1b5Qg2IkjGwYRGAgnrrxZsYJg2X7idrzCJ4Eo32yHms5SSlXAasAFi9enOCCl4ZRQ6p6Wn0wEFHv15sB49GMFcag7fdFLRB+DV4KHByhsRgSZZwXOAxE1Pv1ZsB4NIaUZCMwUwgxHagDLgYuHdkhGeLCFDgMSNSNN2MYk0gpvUKIm4C3UUkAz0gpPx3hYRniwRQ4DEjUjTeT4ozzq9Q+kVK+Abwx0uMwJMhIFziMApIWdePNpDbmKtUwJknlKr9BYkB56sabSV3MVaphTJLKVX6DhKkoHaeYq1TDmCSVq/wGCSPq4xRzlWoYs4yz/NtIjKiPU8xVqsEwNjGiPk4xV6kGw9hESDl89UBCiCbgwCDtrhg4PEj7Gi5SbcypNt6pUsqS4X7TQbbrcFLp+JuxDg16rHHb9rCK+mAihNgkpVw80uNIhFQbc6qNd6yRSsffjHVoSGaslqEajMFgMBiGHyPqBoPBMIZIZVFfNdIDSIJUG3OqjXeskUrH34x1aEh4rCkbUzcYDAZDb1LZUzcYDAZDBEbUDQaDYQyRkqIuhDhLCLFLCLFHCHHnSI8nEiHEFCHE34QQO4UQnwohbg1snyCEeEcIURn4WzjSYw1HCGEVQnwihFgbuD9dCLE+MN7/FEKkjfQYxwOj2b5T0bZTxa6FEAVCiP8SQnweOL6nJHNcU07Uwxa8/jpwHHCJEOK4kR1VL7zAj6SUc4ClwI2BMd4J/FVKORP4a+D+aOJWYGfY/RXAw4HxHgGuHpFRjSNSwL5T0bZTxa4fBd6SUs4GTkCNOfHjKqVMqRtwCvB22P27gLtGelz9jPk1YBmwC5gU2DYJ2DXSYwsbY2nAaE4H1qKWKzwM2KIdd3Mbsu8hpex7tNt2qtg1kAdUEUheCdue8HFNOU+d6AteTx6hsfSLEGIasABYD0yUUtYDBP4eNXIj68UjwE8Af+B+EdAqpfQG7o/q4zyGSBn7ThHbThW7ngE0Ab8LhIqeFkJkk8RxTUVRj2vB69GAECIHeAW4TUrZPtLjiYUQ4ptAo5Ty4/DNUZ46Ko/zGCMljnsq2HaK2bUNWAg8IaVcADhIMoSVil0aU2LBayGEHWX0f5RS/imw+ZAQQgIzUV9aYwL7WwnUSSn/36APFk4FzhVCnA1koC4FHwEKhBC2gFczKo/zGGTU23cftj1JSlkvhJhEArY9hKSSXdcCtVLK9YH7/4US9YSPayp66sEFrwOz1hcDa0ZqMEKI54UQ9UKIdiHEbiHENUIIAfwW2CmlfCjs6eHjvAIVj4wLKeX3h0jQkVLeJaUslVJOQx3P/5FSXgb8Dfh24GkJjdeQNKPKviPpx7avCPw/KmwllexaStkA1AghZgU2nQF8RjLHdaQnCJKcVDgb2A3sBX42wmOZC6QH/p8NNKBm0yWwDdgSuJ2NiudJYD9q8mbCSB/LKJ/nK8DawP8zgA3AHuBl/TnNbci/g1Fj31HG9sU+bPuvQOVotO1UsGvgRGBT4Ni+ChQmc1xNm4BBJHCWfQ+4VUq5OrDtx8APUT+Eu1Fezkwp5Z6I114M3C7D2mwKIX4AfFVKea4Q4lnU5dndgce+CdwHTEOd0b8vpdwmhPgn4Hwp5TmB5+0BNksplwfu1wDnSCm3DM1RMBgMI0kqhl9GHUKI3wghnMDnQD3wRmD7WcDtqJSvmcCZfexmDTBLCDEzbNulwAtR3m8h8AzwPdSZ/ElgjRAiHfhf4DQhhCUQg7OjYosIIWYAOShPwGAwjEGMqA8CUsobgFzgNOBPgCvw0HLgd1LKHVJKB3BPH/twouJllwAExH020eOp1wJPSinXSyl9UsrnAu+5VEq5D+hAXcp9GXgbqBNCzA7cf19K6Y+yT4PBMAYwoj5IBMT1A9Rs+vWBzcfQM+e4vyXPXiAg6igv/dWA2EcyFfiREKJV31AZE8cEHv9fVAzxS4H/30MJ+pcD9w0GwxhluFMax3wA/+qrryY7Oxvg0SuvvJKJEyd+DXgMYPfu3VRUVFBZWVkZ7bUej4dJkyaxZcsWOWvWLB5++GGA7wBcccUVlJaWAvzsuuuuo6ysjJ/97GeZEbt4AXhh1apVvP7665dVVVXx5ptvXr9161b++Mc/sm7dOl5++WWAB4bkw488Ztlsw7jHeOoDoLGxkZdeeonOzk58Ph9vv/02L774IqeffjoAy5cv59lnn+Wzzz7D6XTyi1/8os/92Ww2vv3tb/PjH/+YlpYWli1bFvV51157LStXrmT9+vVIKXE4HPz5z3+mo6MDgC9/+cv87W9/o6uri9LSUk477TTeeustmpubWbBgweAeBIPBMKowoj4AhBA88cQTlJaWUlhYyO23384jjzzCeeedB8DXv/51brvtNk4//XTKy8uDYt8Xl156Ke+++y4XXnghNlv0C6nFixfz1FNPcdNNN1FYWEh5eTnPPvts8PGKigpycnI47bTTAMjLy2PGjBmceuqpWK3WgX9wg8EwahnulMYxH34xjCgm/GIY9xhP3WAwGMYQRtSHAb/fjynyMhgMw0EqNvRKGaSUeL1e3G43brcbm82G1WrFarVisViwWq0IIVDtNAwGg2HgGFEfArSYe73e4DYhBH6/H5/Pp/ozhAm5xWLBZrMFhd5isWCxWIzYGwyGhDGiPohEirn2wv1+P0IILJbe0S7dhMftdlNTU0NmZibFxcUAQZEP9+6N2BsMhr4woj4IaA88UszjIdpzrVZrUOw9Hg9ut7uXZ2/E3mAwRMOIepJo0fV6vfh8PiC2mCci8OH/R3udnnD1er14PJ4ejxmxNxgMRtQTJNyD3rlzJzNmzCAtLW3YhFO/T+T7RYq9lBKfz8ehQ4coLS3FarUG4/bhN4PBMLYwoh4nUkr8fj9erxe/XzU5dDqdvSY9B4IQIunUx2hi7/P5aGlpYcqUKT3CQ+GvieXZG+/eYEhNjKj3g/Z4vV5vUMD1zWKxBAV+NCKljCnS+uQRS+x1KCfcuzdibzCMfoyoxyCamEeGK1JB1GOJcKwwjn6dnvx1u929Hvf7/WRnZwe9eyP2BsPowYh6BOGTn7HEXDPYoj6Q8Es0kg0N9TVJ297eTl1dHbNmzerxWGQYxxRWGQwjgxH1ANFyzPubSExlTz0Z9L60eIe/jymsMhhGB+Ne1GMVDMVDvKI+UgI22KIOKgYf2b63L89eF1ZFE3uTfmkwDD7jVtR1JktzczPNzc2Ul5cnLChjNfzSF36/P+5UyP7EPrKw6i9/+QttbW1ce+21gzpmg2E8Ma5EPVrBkMViweVyJSV+8Yq61+ulu7ubzMzMYfVEh8tTT5RYYl9XV4fdbh/Qvg2G8c64EPVwz1CLsBYWq9WatLdtsVj69Ky9Xi8HDhygvr6e9PR0XC4XVquVrKwscnJyyM7OJicnB7vdPiRiP1SiPlRFS+3t7ZSXlw/Jvg2G8cKYFvVoBUORXqLVag167Ymim3VFosW8oaGB0tJSTj755KAY+nw+HA4HDoeD5uZmqqurg215AdLT08nKyiI7O3vAXutQhV+Gakm8trY2CgoKhmTfBsN4YUyKel8FQ5HYbLZexTfxEhl+CRfzKVOmsHTp0uCVgD5xWK1W8vLyyMvL67Evj8fD/v37cblcHDp0CIfDgdfrJS0trYdXn5WVFbeoDpWox1o7daC0t7cbUTcYBsiYEvV4CoYiGchkp36tx+PhwIEDHDp0iClTpnDKKackHKKw2+1kZWWRlZXF5MmTg5/H7XbjcDjo7OyktrYWh8OB3+8nMzOT7OzsoNhnZmb2es/RGlOPRVtbG/n5+UOyb4NhvDAmRD2RgqFIBhJ+kVLS2NhIdXV1v2Iej7hGhnOEEKSnp5Oens6ECRN6vG93d3dQ7A8fPozT6QQIhm6ys7N7dXEcDBLJfkmUtrY2CgsLh2TfBsN4IaVFXYt5TU0NOTk55OXlJSw4yXiy2jOvra0lPz8/Kc98IAghyMzM7LGgBijB7erqorOzk46ODlpaWuju7qa1tTUo9NqzT7az5FB76kbUDYaBkZKiHrkoRVdXFzabbcgv3cPDLGVlZVRUVNDV1TVqWthaLJagcAPk5ubS2dlJWVkZTqeTzs5Ojhw5Qm1tLS6XC5vN1kPo45mcHcrsF5fLRWZm5pDs22AYL6SUqOu0xMhFKQYy2RkPehKzsbGRsrKyoGfe2NgYdzw+nvj2UBUfWa1WcnNzyc3N7fG4x+MJin1jYyOdnZ3BydlIsdfe+VBlvwzm5zYYxjMpJ+o+n69XJovdbh+S+HEsMdfESmmMZKQaW/UX/7bb7eTn5/e4wtEnzs7OThwOB3V1dTidTnw+HxkZGTidTlpaWpBSkpWVNaheu2kAZjAMnJQS9VgToDabja6urgHtN1wA3W43Bw4ciCnmmrHY0EsIQVpaGhMmTIg6Obt9+/ZgGEpPzupMHO3VJ1M56/F4hixV0mAYT4yJX9FAwy9anHWeeWNjI1OnTu13ArS/itJEGc29X/TkrNVqpaysLCjAenLW4XDQ0dFBQ0NDcJ4hvHI2Ozub9PT0mONpa2vrlbtvMBgSJ6VEPZYgDIao7969myNHjsQl5uGvG62eujh8mMwPPlD/n346MixLZiBExtTDJ2ePOuqo4Hafz4fT6cThcPSYnLVarT28+uzsbNLS0kw1qcEwSKSUqEN0bzZZUXe73ezfv5+2tjaKiooSTk0craJue/llMr7/fbL1PIPNRveqVXi//e0B7zte7z/W5KzX6w22SWhqamL//v14PB4efPBB6urqWL16NcuXL4+6z7feeotbb70Vn8/HNddcw5133tnjcSFEOvB7YBHQDFwkpdwfeOwu4GrAB9wipXxbCDEl8PyjAT+wSkr5aCLHw2AYbaScqEcj0YlSLeZNTU1MnTqVkpISSkpKEp70G42td8Xhw2Rcdx0ivKDK6yXjuutwfOUrg+axJ4tOPY1MP73oootYu3YtaWlpUV/n8/m48cYbeeeddygtLWXJkiWce+65HHfcceFPuxo4IqUsF0JcDKwALhJCHAdcDMwFjgHeFUJUAF7gR1LKzUKIXOBjIcQ7UsrPBvtzGwzDxehIsE6AgfRvcbvd7Nq1i40bN5Kdnc0pp5xCaWkpdrs9KU8/XlHXRVLDgeXPf4ZoFbJeL7bHHhvw/ocqO0UIwfz58/nWt74V9fENGzZQXl7OjBkzSEtL4+KLL+a1116LfNp5wHOB//8LOEOoAZ8HvCSldEkpq4A9wElSynop5WYAKWUHsBOYPPifzmAYPsaEp97fhKXb7aaqqorDhw8zbdo0Zs6c2cMrT7b9bjyifuTIESorK/F4PEgpSU9P79WgS49loJ667eWXybj55qiPCSD9scfw3nTTgLz1ocon72+itK6ujilTpgTvl5aWsn79+sinTQZqAKSUXiFEG1AU2P5R2PNqiRBvIcQ0YAHQa6cGQyoxJkQ9FlrMm5ubmTp1ai8x1+iWuInSl6i3t7dTWVmJxWJhzpw5wbCC2+0OlvE3NzcH0wJ1FagQArfbHTMMEQtx+DAZV19Nn3603Y44cGDEQzDRaGtro6ysLObj0U4mUa4aon182cd2vZ8c4BXgNillexzDNRhGLSkn6vFc/rtcLvbv309zc3NUzzySZCdao4m6w+Fgz549uN1uZs6cSUFBQbDbYniDrqKiouBrdKbIwYMH6ezs5NNPP8Xj8QTb7upbX8U+lrVr+x+ww4Fl61b8ixYl/FlhaJt5tba2Mn/+/JiPl5aWUlNTE7xfW1vLMcccE/m0WmAKUCuEsAH5QEvY9uDugIMAQgg7StD/KKX808A/icEwsqScqMdCCEFXVxcHDhygpaWFadOmUVFREXemRrLhF+1Bdnd3s2fPHhwOB+Xl5T1EO573z83NZcKECaSlpTF9+vTgiaCzs5POzs5eXn242KelpWHZurXf9xFAxk9+guPcc5Py1odS1Pvrpb5kyRIqKyupqqpi8uTJvPTSS7zwwguRT1sDXAGsA74N/I+UUgoh1gAvCCEeQk2UzgQ2BOLtvwV2SikfGoKPZTAMOykn6tFE2uVy4XK52Lx5M9OnT2fWrFkJTehZLJakPHUhBD6fj88//5wjR45w7LHHUlJSMiiTibG8er/fH2y529zczIEDB/B4PEwHpsezY7cby9at+M44I+ExDXWHxr5E3Waz8dhjj/G1r30Nn8/HVVddxdy5c/n5z3/O4sWLOffcc0EJ9B+EEHtQHvrFAFLKT4UQq4HPUBkvN0opfUKILwLfBbYLIbYE3uqnUso3huRDGgzDQMqJejgul4uqqipaWlpIS0vjhBNOCMamE8Fms+FyuRJ6jdfrpaqqCqfTSV5eXsInkmSxWCy98r+llNieeirufSQ71TnUnnp/bXfPPvtszj777B7b7r333uD/Uspu4MJor5VS/ivwrxHbPiB6vN1gSFlSTtSFELhcLvbt28eRI0eCnvmOHTuSXuwikYUyfD4f1dXVHDx4kClTppCdnR0ttpsUyWa/WJqbyXz55bjVaXddHezZQ3Z2Nrm5uXE35hpJT91gMMRHyol6W1sbW7ZsYfr06cyePTvoHQ90rdH+RN3v91NXV0d1dTXHHHNMcP3R2trapN5zMLFs3QoJzAmUS0nLhAl0dnZSXV2Nw+EACPZq0bf09PQerxtKT93pdJKVlTUk+zYYxhMpJ+p5eXmccsopvUIdAxF1m80WU9SllNTX17N//35KSko46aST+l1IYrhJ1LcXS5f26sLo9/t7LKRRU1OD2+3GbrcHRX6oRF1fnYyWxUYMhlQm5UTdYrFEjV0nWxUK0cMvUkqamprYu3cvBQUFLFq0qJfnmgjxrlGaTPhFnnAC2GwQz+cXAqJk5lgslqB4hxOZgeNwOGhra+vl1Se7PF7PoZnwtsEwUFJO1GNhs9mSXigjMqWxpaWFyspKsrOzOfHEE0f9EmvWv/0NfL6YVTY9kJL0H/6Q7t//Pq59h/dW14tkRC6PF82r13n18cTgvV7vkMXqDYbxRsqJel/tdxPNYNFYrVa8Xi9tbW1UVlZis9mYO3duL681FoPZtzxRxOHDZNx4IyJOD18AtldfRezahZw1K6H30uuTxuPV19TU4HQ6gysk9eXVt7e3m17qBsMgkXKiDoPbfhfUwtXt7e3s2bOHioqKhARGFyD1J+pDtUapOHAAkrhCsW7ahDdBUe9vfdJoKybpRTQ6OjpobW3tseh1bm4u2dnZ7Ny5s89j3l/LXZfLxeWXX87q1av3EEfL3cD2Z4BvAo1SynkJHQiDYRSTkqIejWRE3el0snfvXpxOJ3a7nUVJlM9HLoXXFx7PYTo6/o7TuQOLJYvCwnPIzExMWCORHk/0roz9IFpaEn6N3+9PeMm58EU0wtHroHZ2drJ69WrWr1/P8uXLWb16dY/nxdNy97e//S2FhYXE23JXSukDngUeQ/VTNxjGDGNK1OONqbtcLvbu3UtbWxvl5eUUFxezbt26pN43nk6NHo+HXbvuxeF4sMf2nTt/Q13ddUyceBOLFuWSTPKHLbC6UTQifX4R9jf9F7/Ae8klCbUL8Pl8A5osDsdut1NYWEhhYSGXXHIJJSUlPPLII72eF95yFwi23A0X9ddee4177rlH3/0v4LHIlrtAVaDS9CRgnZTy74HOjAbDmCIlRT1amCKe7BePxxNswTtjxgzmzJkz4Fh4X6Lu9/upqtpCQ8MPsFh6dnRds+ZaHn30cfx+W2A/fh59NJsvfKE5ofe39u4p3gM3qml+ry86iXYBQ5XSqAuPon0X8bTcDX9Ooi13DYaxRkqKejT6Cr94vV6qq6upr69n6tSpLF26dNDEKZqoSylpaGhg376nsFh+hcXS8wrixRd/wKpVDxKeq+L3W7n55gm8/PKnZGXVBScV+4phi127sG3d2mfGS19+daLJk0NVUdpXNWk8LXdjzEP023LXYBiLpKSoR/Pooomr3++npqaG2tpaJk+eHKwCHUwi3/fIkSPs3r2b7GwvVuuvkLKnoP/nf/YW9BCC3/9+KitWdHHw4EEcDgd+v5+sLDeZmW3k5VVQUDAtWPxk3bRpQGNvyc4mkcL8ofLUW1tbmTw5ugMdT8td/ZzS0lLibblrMIxVxkwJX7jQSympq6tj3bp1eL1eTj75ZKZNm9a31xuY8EwULeqdnZ1s3ryZqqoq5s6dS0nJrl6C/tpr17JyZSxBBxCsXVtKRsYUystLmDNHMHXqdtzu82lru5bq6i+wbdsjbNiwgW3btlGTkRFzXPG4o9V9xOOjMVSeekdHR8zsl/CWu263m5deekl3ZAxy7rnn8txzehW7UMtdVCvei4UQ6UKI6QRa7g76BzAYRhFjxlOH8LDHPoqKiliyZEncKwjpAqREPVEpJXv37sXlclFRUUFuro9duy6io+OvPZ7X2lrMf/zHr+mvPEhKK7/85Wa+852zEMKK39/Z43EhHuDEE6/E58vB9sorfe6rv9mC5j/9CW6/vZ9nhRjKmHqsDo3xtNy9+uqr+e53v0u8LXcBhBAvAl8BioUQtcC/SCl/O+gfzmAYZlJS1CORUgYXkWhubmbhwoVk9OHFRkMXIMWbsufz+di/fz9NTU2UlZVx4okncujQU2zadEvU569Zcx0+XzyZI4Knnz6NvLzLOOecp3s/Kmy4XAfI6Soj5ze/6cPn7+9d4KwdO3j217/mkluijzmSkYipQ/8tdzMyMnj55ZcByiNfG63lbmD7JcmP2GAYvaR8+KW1tZVNmzZx8OBBcnJyqKioSFjQIf72u1JKamtr+eijj4A2Jk48Qm6un9rah6iqii6Ora3FPP/83cTfulvw8MMraW3tnW7o97uBTP76zDN4k6yg1diBl+++m/fffz+u54+Ep24wGBIjJT11IQQdHR1UVlYCMHv2bHJzc9m8eTNerzepLorxLGl3+PBhKisrKSws5NhjD3DgwI34/VY6OrpQBYvRefjh/8DjSexEI6WFv/zlUpYv/3Vgix3wIISFysqv0vFSN4PRK/IrwDe+8Q3OO+88br/99uACHDk5OWRkZPQIdQ2lp25E3WAYHFJS1Nvb2/n888+DCztrBtIqQIdfotHR0cGuXbuw2+2ccMIJ2O0ONm++Eb+/q9/9Hjgwi7///SKSWWBn5coHKSo6zJlnvhycdJWyC+cBuKRq4Ev2COAHwCOoAp5p06YFwxz19fV0d3cHy/lzcnLweDxJdZHsD4fDEXefHYPB0DcpKer5+fksWbKk1/aBinpk+KW7u5vKykq6urqYNWsW+fn5SClpb9+B8pz7F/X33z8vqfGAQEob//Zvz7Jo0V8oKDgMwLvvwv+sgC9LBsVT9wDTgMPAo48+yqOPPsp3v/tdHn/8caBnky6Px8PmzZsRQgRXTdKCn6wHL6VESml6qRsMg0RKinos7Hb7oLTf9Xq97Nu3j8OHD1NeXk5JSQmgwg9+vx+7vbRXumIsqqsrkhqPxuez8cEHX6a8/BUyM+FXv4LpXohcIyiutrtRsAP7I7b94Q9/YObMmdxyyy3YbLZgk65Dhw6xZMkSfD5fUOjr6+vp7OwM5NNnBUU+Nzc37swjML3UDYbBYkyJ+kA9dY/HQ3V1NdXV1ZSVlbF06dJgSwK/3x/stJiePpHp059g//7rgd5ph5pHHnmUd965agCfSPHgg+eTkfEKLhdI2bOaJll0EOUplJceyb333stFF11EcXExPp8veAx8Ph9CCPLz88nPzw8+P3zlpJaWFqqrq3G73WRkZPQQ+sg4vd/vN4JuMAwiKXnN21dP9WREXUqJw+GgsrISl8vF0qVLKSsrQwiBz+cLhmXCV10qLr6QE07YydSpDyJEdq99bt++lNdeuxnlPw9EtASwnO7uWehw9g2xPkfErb+9CuB6IFpLL6/Xyx//+EfS0tJwu93s2LGD4uLioLB7PB7cbjderzd4fHJycjj66KOZOXMmCxYs4KSTTmLWrFnk5ubS2dlJZWUlGzduZPPmzVRWVlJfX8+OHTvIzc2NOc633nqLWbNmUV5ezv3339/rcZfLxUUXXUR5eTlCiPXhTbqEEHcJIfYIIXYJIb4Wtv2swLY9Qog7e+3UYEhhxFBMfPXBoL2Z2+3uNWlXX19PV1dXsKNfPLS2trJ7926klEyYMIGZM2cGvVIdjhFCxDyReDxNbN06p8ek6bvvXsz99z+Lz5fGwKczQR02H/AHZnEvO9k/KHvVe/4lcHeUxzIyMnjzzTcBlWGkqz71VUv430j0CTBarNzj8dDR0UFnZycPPvgga9as4Ytf/CKvvvpqj+f5fD4qKip6tN198cUXe3Ro/M1vfsO2bdtYuXIlQohLgH+UUuq2uy+iujIeA7wL6FjYbmAZqo3ARuASKeVn8R81g2H0Mm7DL06nk927d+Pz+ZgzZw4OhwOHwxGMm0PfYq6x20soKfkVDQ0/QggbBw5MZcWK38VZaBQvAvVVXcnJCOCfBnXPP0JlwESGYYQQNDY2cvbZZ/cQZ/1/+ORopNCHh2tAXQ1ZrVaEEFit1mCc/qqrrsLr9fLMM8/0GtsQtd0F2COl3Bf4jC8FnmtE3TAmSFlRT3b1I7fbzd69e2ltbaWiooKioiKklHR1ddHY2IjdbicvL4+cnJx+BV2HbCyWRRx33Fb++7/93HTTLHy+eP3oRKc3Bes5OYHnxz+KORkZvN/d3WN7V1cXhw4diiszJZbQ67+RXr2O0zc3N5OWlha1kncI2+7WRGwf/INqMIwQKSvq0egr+8Xn81FdXc3BgweZNm0as2fPBkKCU1BQwLHHHktHRwdVVVU4nU6sVit5eXnBW1ZWFv66OlYAABJ3SURBVEIIPB4P+/bto62tjZkzZ1JYWEhTE/zgB5kJCfrixV1s2ZKO12uhp7jHFvsFbMGPmgwZrBBMBvC7K6/kqfx8VqxY0eOxO+64g3PPPTeYAZQIWugjTwpa4Hfs2MHdd9/NF77whaivH6K2u9HOUKYdr2HMMKZEPZqnLqWkvr6eqqoqJk2aFGy/K6UMeos6JFBUVERRUVHwtR6Ph/b2dtrb22lsbMTpdOL3+/F4PEycOJG5c+eSlaWSC6urLURP1dZ60VO0ly/38rvfQVOTi2eesfHv/25HSnC5QAgtVj11qZgmnuc7DHZNpwCmr1rFcf/+7+Tk5NDZGcrmsdvtVFdXJyXqsfD5fDz88MO8+eabPPfccyxcuDDq84aw7a5px2sYs6Rk9gtEz4CJFPWWlhbWr19Pa2srS5YsYcaMGVgsFnw+H16vNyjoscIsdrudoqIipk+fTmlpKRaLheLiYubMmUNGRgZ79uzho48+4uOPP8br3UN3dzSHL3rf9FdftfHuu+rw33GHl7Vru1CRCYGU0TJmJP+Pu7CQeHvgeBB+PwsDVyHhuN1uCgsLk2pLHI3t27fzta+pRJT3338/pqDDkLXd3QjMFEJMF0KkoTo6rhmUD2cwjALGlKeuq0I7OzvZvXs3QgiOP/54srOzg555IpOgoOLmu3fvxmazMX/+fDIzM3s9x+VyUVXVicUS/xrQbjdcemk6fj9cfrmXZ5+10VfdVDFNXMMzgxZyiUbFn//Mk08+yfXXX4/dbsftdnPffffR1dXFxo0bAYJVpHreId5KUrfbzYMPPsg777zDypUrOfHEE/t9TSJtd8vLywF+SHxtd28C3gaswDNSyk8TPVYGw2glZVMatbcdjsvl4oMPPgh2awysMB93emI4Om7e3t7eq8dMND7+2MI//EM63d3JyG5/E6aSr4m/8KY8a0hFXQLeK66g4ZZbOLxpE8WLF1MUmHsAgifM9vZ2Ojo66OjoAAgWFmmhj5z03LZtG7fccgvf/OY3ufPOOxOqNE0QU8VkGPeMCVH3er3s37+fQ4cO4fV6+dKXvqTeLKISNB4x9/v91NbWUldXx7Rp0zj66KP7fV1TE/z97xYuvzyd/ic8o8fY+9Kjo62NPCpv5kL/6iFXLQmQnq5uHg/uJ57Ad+GFMZ+vV30KF3q/3092djZbtmxh3bp1fPbZZzz99NPMnz9/iEdvRN1gSFlR9/v9uN1uamtrOXDgAKWlpZSVlfHRRx+xdOnShMUcVGvdPXv2UFJS0u/yd5rVq63ccEMaFgs4HNCfrlit6hD0lyVTTBPTqOIk62YeET/A5u3uV7EiTw0SkGrWNe661l77SEuja906CPPY+8Pv97N582buu+8+nE4nEydO5JV+VmkaJIyoG8Y9KSvq7e3tbN68mQkTJjBjxgzsdjtSSjZs2EBxcTH5+fnk5ubGtZKRjsGnpaVRXl4e9yIbTU0wZ04mXV19h040aWmwYoWTn/wkE48ndmbdd6wv8KTvGixWSPf1L+a99pCeDlLi+tWv2L9gAZ07dlCxZQu5Tz8NOgzVx+sjTwzYbLiffrpPj13jcrlYsWIF77//PitXruT4449PYPQDxoi6YdyTsqLudrtxOp1kZmb2iJt3dXVx5MiRYDhACBGM9+bn55OdnR303N1uN/v27aOjo4OKiooeDari4eOPLXzzm+m0t/cVSpHccIOHf/gHH0cdVcfmzS3cddciOjp6XwWkp0v+8sdqvnjhsVhknDOuEcicHNz33ENnTg4H6+rIKS2lZNkyLBMnYn36adJ+/GNwuxNWP5mRQdfnn0MfqY2bN2/mtttu44ILLuD2229ParGSAWJE3TDuSVlRl1Licrn6nQT1+Xx0dHTQ1tZGe3s7DocjWK7e1dXF1KlTg+mKiRLNU7fZJCrUH9qWkeHnxRfXMX16Drm5M5g3LyfCu5dkZMDKlW4uW30+1jfeSFqd/FarCreEpyDabLgfeoi0O+5AdIV61ESL7sdCZmfjevNN/IsW9Xqsu7ub+++/n//7v//jySefZO7cuUmOfsAYUTeMe1I2T33Lli288MIL7NmzJ7jIQrTYudVqpaCggKlTpzJv3jyOPfZYvF4vGRkZTJo0iebmZtavX88nn3wS7KHudrvjGkNJCTzxhJvMTEleniQzU/LTn3qIbDpotfrJyjqOmTNncvTR1h6vyciQ/PznHj7/vIsLj/80KUEP78wofD4sfn8whi4A4fWS9oMfqPh6OJmZ+M46K74zrc+Hv6ys1+ZNmzaxbNkyiouLee+990ZS0A0GAynsqe/evZuXXnqJjRs3UlVVxeTJk1m8eDGLFi1i8eLFFBUV9RB5vaZpeno6xx57bI+4uZSS7u5u2tvbgx691+slOzs7GLbJzc2NOXHa1KQqSsvKlHcc6b1nZkp27uzqEbkIf43ebn3+edK+970hcTdj5dxgtYLP12eMHbsd91NP9Yipd3d382//9m+sX7+eJ598kjlz5gzBqBPGeOqGcU/Kino4fr+f6upqPvroIzZs2MDGjRtpa2tj9uzZzJo1i+bmZpYvX96jfWy/Aw30WNci39HRgZSyV3w+MmzT2trKE08c4YEH5pCWJvB4lDd/4YVxxMg//5zMRYvizlKB+FUs3tZhQaHPzAS3G8/NN+O9+eYesfSNGzfywx/+kEsuuYTbbrstrsnoYcKIumHcMyZEPRoej4eNGzdy6aWXMn/+fGpra7FarSxYsICFCxeyZMkSKioqElpbU8fndT+Yzs7OYNOvrKwsmpub8fv9zJo1C6czu5cnHg/2H/0I28qVwfvRstyHorxJ48vOpvKXv6TDZsMzeTIZU6aQl5dHRkYGdrudFStW8PHHH/Pkk08Gm6KNIoyoG8Y9Y1bUNX6/H4vFgpSSzs5OPv74Yz766CM2btzI7t27KS4uDoZtTjrpJCZOnJjQ8moul4u9e/fS1NRERkYGfr+f9PR08vPzgx59XBWUTU1YqqtV3Lq5GeumTXDgAGkrVgR7D8RqOxj5WKySp7g89cxMunbuhJISvF5v8AS2du1aHnzwQXJycnjxxRc56aST+t/Z8GNE3TDuGfOi3he6g+OGDRuCQt/Y2Eh5eXkwNr9gwYKYvdVbWlqorKykuLi4R7FSZHze7XYH4/P6Fh6ysK5eTdoNN4Dd3ruKs6kJ6yuvkPbP/4xwOnt/howMPP/0T9ifeKLPuLhv4UKsmzf3HTtPT8f95JM9YudOp5P77ruPTz75hKeeeoqSkhLsdjs5OTn9H+Dhx4i6YdwzrkU9Gj6fj127drF+/fpgVozH42H+/PlBoc/JyaG1tRWr1UpFRUWw/W4spJQ4nc6gyLe3t+P3+8nNzaXQ6+XY009HhC1QEe4tA7hqa8mbPx+ryxV6DkBGBu6VK5EzZpD+jW8gAr1Yer0/QE4OdHbGFvX0dLr+7/96VI6uW7eOH//4x1x++eXcfPPNCYWqRggj6oZxjxH1OHA6nXzyySesX7+el156ib1793LaaacxdepUFi9ezJIlSxLOdff7/XR0dOD+8EMmX3klNtVjAABfTg7tr7yC/QtfoLa2loMHDzL/s88o+slPlDfvduP5yU/wXnWVEv6mJjLnzOmdg56dDQ5Hn6GZ8JOD9tAdDgf33nsvO3bsYNWqVcycOTOJozYiGFE3jHuMqCfI9u3bKS8vx+Fw9Ajb1NTUUFZWxpIlS1i0aBGLFi2ioKCg//h8FEH2Z2SwYfVqmoC0tDQmTpxIQUEBBR4PGQ0NKu4eMftqffll0q6/PhTCWbEC0tNJu/32Hh68zMhQMfqMDPB4epwcpJR8+OGH3HHHHVx11VXccMMNw+Kd+3w+Fi9ezOTJk1m7di1VVVVcfPHFtLS0sHDhQv7whz/E29nRiLph3GNEfZDw+/3s27cvGLbZtGkTDoeD4447jsWLF7N48WLmz59PenrvBakjBfnAz3/OgVNOYc6cOdjt9h7xeZfLRVZWVo/8+WA5fvhkaywPPjOTrg8+wOJw9Dg5OBwO7rnnHj7//HNWrVrFscceOyzHDeChhx5i06ZNwQnZ5cuXc/7553PxxRfz/e9/nxNOOIHrr78+nl0ZUTeMe4yoDyFut5tt27YFhX779u2kpaWxYMGCoNCXl5ersE1TE+3btlHp8TBx3jwmT54c1cvXi2SHx+d9Ph85OTnBjJvc3NzQQtCRHnxEK10pJR988AF33HEH1157Lddff31SLROSpba2liuuuIKf/exnPPTQQ7z++uuUlJTQ0NCAzWZj3bp13HPPPbz99tvx7M6IumHcY0R9GJFS0tbWxqZNm1i/fj0bNmxg7969TJgwge7ubr73ve/xpS99iWOOOSahtMrwnuZtbW10dHRgsViCmTb5bjc5hw8jp07tEbbp7OzkX/7lX6isrOSpp55i+vTpQ/Gx++Tb3/42d911Fx0dHTzwwAM8++yzLF26lD179gBQU1PD17/+dXbs2BHP7oyoG8Y9o6YUcDwghKCgoIAzzzyTM888E1BCf+ONN5Kfn8/OnTv5/e9/T0tLCxUVFUFv/sQTTyQrKyum0IcLeGlpKaAWDtGNzPa2t+PweLAfOEB+ayvNzc3U1dXxwAMPcP311/P4448Pq3euWbt2LUcddRSLFi3ivffeA/SC2z1J5ARnMIx3jKiPMEIIfvOb3/TY5vV6+eyzz1i/fj2rV6/mrrvuAuCEE04ICv2sWbP6LM+32WwUFhZSWFgY3OZ2u2loaOB3v/sd7733Hjk5OVx66aUjIugAH374IWvWrOGNN94I5vbfdttttLa24vV6sdls1NbWcswxx4zI+AyGVMSEX1IA3Yfm448/ZsOGDWzYsIFdu3ZRWFgYzJ1fsmRJn2EbKSXvvfceP/3pT7nxxhu55pprEloVaqh57733eOCBB1i7di0XXnghF1xwQXCidP78+dxwww3x7GZ0fBiDYQQxop6iSCk5dOgQGzZsCMbn6+vrmTFjBosWLWLJkiUsXLiQ3NxcOjo6uPvuu6mtrWXVqlWURWmhO9KEi/q+ffuCKY0LFizg+eefj5o1FAUj6oZxjxH1MYTf72f37t1Bkd+8eTPt7e20trbyi1/8gquuumpIQy01NTVcfvnlNDQ0YLFYuO6667j11ltpaWnhoosuYv/+/UybNo3Vq1f3CAsNIkbUDeOelBf1t956i1tvvRWfz8c111zDnXfeOdhvkdJ0dXVRWVnJ/Pnzh/y96uvrqa+vZ+HChXR0dLBo0SJeffVVnn32WSZMmMCdd97J/fffz5EjR1ixYsVQDMGIumHck9Ki7vP5qKio4J133qG0tJQlS5bw4osvctxxxw3m2xiS5LzzzuOmm27ipptu4r333mPSpEnU19fzla98hV27dg3FWxpRN4x7UnY5O4ANGzZQXl7OjBkzSEtL4+KLL+a1114b6WEZgP379/PJJ59w8sknc+jQISZNmgTApEmTaGxsHOHRGQxjl5QW9bq6OqZMmRK8X1paSl1d3QiOyACqqOmCCy7gkUceiXulKYPBMDiktKgPRqFKTU0NX/3qV5kzZw5z587l0UcfBVSv9GXLljFz5kyWLVvGkSNHBmXMYx2Px8MFF1zAZZddxvnnnw/AxIkTqa+vB1Tc/aijjhrJIRoMY5qUFvXS0lJqamqC95MpVLHZbDz44IPs3LmTjz76iMcff5zPPvuM+++/nzPOOIPKykrOOOMM7r///sEe/phDSsnVV1/NnDlz+OEPfxjcfu655/Lcc88B8Nxzz3HeeeeN1BANhjHPcE+UDipCCBuwGzgDqAM2ApdKKT8dwD5fAx4L3L4ipawXQkwC3pNSzhqEYY9ZhBBfBN4HtgP+wOafAuuB1UAZUA1cKKVsGZFBGgxjnJQWdQAhxNnAI4AVeEZK+a8D2Nc04O/APKBaSlkQ9tgRKeWQJFcbDAbDYJHyoj5YCCFygP8F/lVK+SchRKsR9RBCiLOAR1Enz6ellCYeZTCMQlI6pj5YCCHswCvAH6WUfwpsPhQIuxD4m3AenhDCKoT4RAixNnB/uhBivRCiUgjxn0KIuJbzGWmEEFbgceDrwHHAJUIIUwxgMIxCxr2oC5Uu81tgp5TyobCH1gBXBP6/AkgmAf5WYGfY/RXAw1LKmcAR4Ook9jkSnATskVLuk1K6gZcAM9tpMIxCxr2oA6cC3wVOF0JsCdzOBu4HlgkhKoFlgftxI4QoBb4BPB24L4DTgf8KPOU54FuD8xGGnMlATdj92sA2g8Ewyhj3/dSllB8Qu7z8jAHs+hHgJ0Bu4H4R0Cql9Abup5IwRjs+ZjLGYBiF/H9gwMtABK4JNAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 画图\n",
    "plotShow()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考博客¶\n",
    "主要参考统计学习方法这本书，书籍地址：http://www.dgt-factory.com/uploads/2018/07/0725/%E7%BB%9F%E8%AE%A1%E5%AD%A6%E4%B9%A0%E6%96%B9%E6%B3%95.pdf。\n",
    "\n",
    "\n",
    "[统计学习方法-代码解读](https://github.com/fengdu78/lihang-code)\n",
    "\n",
    "[EM算法 - 期望极大算法](https://blog.csdn.net/randompeople/article/details/93711747)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
